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Preface 


The  stratified  concept  being  presented  here  was  originally 
theorized  by  Col.  Ronald  Bagley  at  the  Air  Force  Institute  of 
Technology.  Captain  A1  Lesmerises  developed  the  details  of  the 
theory  in  his  Masters  Thesis  con^leted  in  March  1992  [2].  This 
work  is  a  continuation  of  that  development. 
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Abstract 


This  study  continues  the  investigation  of  the  stratified 
plate  model  for  the  analysis  of  composites  plates.  The  strata 
theory  models  a  composite  plate  as  a  stack  of  isotropic  matrix  and 
orthotropic,  fiber  dominated  layers.  Energy  methods  are  used  to 
derive  a  set  of  partial  differential  equations  and  boundary 
conditions  as  a  function  of  the  displacement.  An  assumed 
displacement  field,  of  a  polynomial  in  z  and  a  Sine/Cosine 
function  in  x  and  y,  generates  a  solution  to  a  simply  supported 
con^osite  plate  with  cross  plies.  A  cylindrical  bending  exanple 
is  also  solved  with  a  polynomial  in  z  and  Sine/Cosine  function  in 
X  and  y,  displacement.  A  Ritz  polynomial  x  and  y  displacement 
function  is  used  to  generate  a  solution  to  a  conqposite  plate  with 
angle  plies.  The  strata  model  is  expected  to  be  useful  for 
predicting  stresses  in  composites  ^ere  the  matrix  material  has 
become  viscoelastic. 


A  STRATIFIED  PLATE  MODEL 


FOR 

COMPOSITE  ANALYSIS 


I .  Introduction 


Composites  Background 

Composites  currently  in  use  today  generally  consist  of  very 
strong  and  stiff  continuous  fibers  surrounded  and  held  together  by 
a  softer,  more  ductile,  matrix  material.  The  fiber  is  the  main 
load  carrying  material  and  the  purpose  of  the  matrix  is  hold  the 
fibers  together,  transfer  load  from  fiber  to  fiber,  and  protect 
the  fibers  from  the  environment.  A  conqposite  ply  or  lamina 
consists  of  a  single  layer  of  aligned  fibers  held  together  in  a 
flat  sheet  by  the  matrix.  Most  con^site  structures  are 
manufactured  with  con^site  laminates  composed  of  many  plys  or 
laminae,  layered  together  with  different  fiber  orientation. 

Because  of  the  unique  and  almost  infinitely  varicible  arrangement 
of  stiff  fibers  and  soft  matrix  the  mechanics  of  ccnoposite 
materials  is  difficult  to  model  mathematically  and,  except  for  a 
few  single  cases,  exact  closed  form  solutions  to  the  differential 
equations  of  motion  have  not  been  found. 


1 


Current  Composite  Analysis  Methods 


The  axis  system  used  for  the  plates  in  this  work  is  as 
follows : 


a 


Figure  1.  Coordinate  System  of  Plate 

The  X,  y  and  z  direction  are  the  primary  coordinates.  The  1 
direction  is  the  off  primary  axis,  fiber  orientation  and  the  2 
direction  is  transverse  to  the  fiber  direction. 

The  state  of  stress  in  a  plate  is  related  to  the  strain  by 
the  generalized,  Hooke's  law  stiffness  matrix  [S]  as  follows: 


1 

r  1 

1 02 1 

1  ^2  1 

1 03 1 

1  £3  1 

K  ^ 

I  I 

|Y23  I 

|T^13  I 

|Yi3  I 
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The  6x6  stiffness  matrix  [S],  coefficients  are  a  function  of  the 
twelve  3-D  engineering  constants  tensile  moduli  (Eu,  E22,  and  E33), 
shear  moduli  (G12,  G13,  and  G23)  and  Poisson's  ratios  (V12,  V13,  V23, 
V21  f  V31 ,  and  V32 )  •  The  subscripting  on  the  stress ,  strain  and 
tensile  and  shear  moduli  has  the  first  subscript  being  the  load 
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face  and  the  second  being  the  load  direction.  The  subscripting  on 
the  Poisson  ratios  has  the  first  subscript  being  the  load 
direction  and  the  second  subscript  being  the  deformation  direction 
(ie.  v21=ei/e2  with  the  load  in  the  2  direction).  Due  to  symmetry 

of  the  stiffness  matrix,  the  Poisson  ratios  are  related  by  the 
relationships,  V2i=(E22/Eii)*Vi2  which  reduce  the  number  of  required 

constants  from  twelve  to  nine  [9:10]. 

The  classic  laminated  plate  method  of  composite  analysis 
involves  idealizing  the  con^site  ply  consisting  of  fibers  and 
matrix  as  a  homogeneous  transversely  isotropic  plate.  That  is,  a 
plate  with  three  planes  of  elastic  symmetry  (orthotropic)  in  which 
one  of  the  planes  is  isotropic. 


Figure  2.  Classic  Laminated  Plate  Idealization 
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Five  engineering  constants  (En,  £22/  V23,  G12)  are 

required  to  derive  the  generalized  Hooke's  law,  6x6,  stiffness 
matrix  coefficients  of  a  transversely  isotropic  plate  [9:11].  The 
engineering  constants  are  derived  from  volume  fraction  averaging 
of  the  homogeneous  properties  of  the  fiber  euid  the  matrix.  The 
fiber  and  matrix  moduli  and  Poison  ratios  combine,  in  proportion 
to  their  volume  fraction,  in  series  or  in  parallel  to  generate  the 
averaged  con^xasite  engineering  constant  [1:  61-86].  The  stiffness 
matrix  [S]  coefficients  are  derived  from  the  5  engineering 
constants  [9:  11-12].  Most  classic  laminated  plate  models  reduce 
the  Hooke's,  6x6  stiffness  matrix  to  a  3x3  matrix  by  assuming  each 
ply  is  in  a  state  of  plane  stress.  This  reduces  the  number  of 
required  con^site  engineering  constants  to  four  (En,  £22^  and 
G12).  The  entire  laminate  stiffness  is  the  summation  of  the  ply 
stiffnesses.  With  the  stiffness  known,  the  strain  on  each  ply  can 
be  calculated  as  a  function  of  the  strain  on  the  mid  plane,  the 
curvatures,  and  the  distance  from  the  mid  plane.  The  strain  on 
the  mid  plane  and  the  curvatures  are  calculated  from  the  inverse 
of  the  composite  stiffness  and  the  applied  load  and  moment  vector 
[1:  190-194].  This  method  has  proven  useful  and  accurate  in 
predicting  the  performance  of  coo^site  laminates  at  low  to 
moderate  operating  temperatures  [ 1 ] . 

One  dra\d>ack  of  this  approach  is  that  the  area  of  each  ply 
where  the  softer  matrix  is  dominant  is  assigned  the  averaged 
higher  composite  stiffness.  Epoxy  and  thermoplastic  matrix  moduli 
are  usually  on  the  order  of  2-4  GPA  \diile  graphite,  fiber  moduli 
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are  on  the  order  of  200-550  GPA  [1:  31-33;  23].  The  difference  in 
stiffness  is  even  more  apparent  v^en  the  leuninate  is  subjected  to 
high  temperature.  The  matrix  materials  generally  have  glass 
transition  ten^ratures  in  the  range  of  50-200  degrees  Celsius  [1: 
28].  Even  high  heat  resistant,  thermoplastic,  matrix  materials 
have  a  maximum  continuous  use  temperature  in  the  range  of  150-300 
degrees  Celsius  [3:  156].  Glass  fiber  materials,  on  the  other 
hand,  have  a  softening  point  on  the  order  of  771-970  degrees 
Celsius  [7:  2. 4-2. 5]  and  properly  treated  carbon  fibers  have 
operating  ten^ratures  of  over  1500  degrees  Celsius  [8:  62-68], 

The  large  difference  in  stiffness,  undoubtedly  plays  a  role  in  the 
actual  behavior  of  the  material.  Under  high  ten^erature  loading, 
due  to  the  softening  of  the  matrix  material,  matrix  dominant  areas 
must  experience  higher  displacements  while  fiber  dominant  areas 
will  experience  lower  displacements  but  higher  stresses.  The 
averaging  of  properties  within  a  ply  can  not  address  the 
differences  in  displacement,  strain  and  stress.  The  fact  that 
most  models  assume  each  ply  is  in  a  state  of  plane  stress  also 
eliminates  the  ability  to  predict  through-the-thickness 
displacements,  strains  and  stresses. 

Debonding  or  delamination  of  plys  is  a  ccxmoon  type  of  failure 
in  con^sites  and  can  be  attributed  to  through-the-thickness 
stresses  and  strain.  The  classic  Lcuninated  Plate  Theory  assumes 
that  each  plate  is  in  a  state  of  plane  stress  and  does  not 
generate  any  through-the-thickness  stress.  This  is  another 
drawback  of  classic  laminated  plate  method. 
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strata  Concept 


As  an  alternative  to  modeling  each  ply  as  a  homogeneous 
transversely  isotropic  plate,  the  strata  theory  breaks  each  ply 
into  matrix  and  fiber  dominated  layers,  each  of  vdiich  will  be 
referred  to  as  a  stratum.  This  idealization  requires  the 
assun^tion  that  the  fibers  are  evenly  disperesed  and  alligned  in 
each  ply. 


Matrix  Stratum 
Fiber  Stratum 
■Matrix  Stratum 


In  order  to  clearly  distinguish  between  the  matrix  and  the 
fiber  strata,  each  fiber  is  idealized  as  a  square  cross  section 
which  has  equal  area  to  the  actual  fiber.  Each  matrix  stratum  is 
entirely  matrix  material  and  has  the  properties  of  a  homogeneous 
isotropic  plate.  Each  fiber  stratum  is  a  combination  of  matrix 
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and  fiber,  dominated  by  the  stiffer  fiber.  As  will  be  shown  in 
the  following  sections,  this  idealization  allows  the  matrix  strata 
to  experience  different  displacements  than  the  fiber  strata.  This 
model  is  expected  to  more  accurately  predict  the  displacements, 
strains  and  stresses  of  composite  plates  with  dissimilar  matrix 
and  fiber  constituents. 
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II.  The  Strata  Modal 


Strata  Laminate  Idealization 

Each  lamina  in  a  laminate  of  con^)osite  plies  is  idealized  as 
a  stack  of  matrix  and  fiber  strata  as  shown  below: 


Figure  4.  Strata  Laminate  Idealization 


The  lower  matrix  strata  of  one  ply  ccanbine  with  the  upper  matrix 
strata  of  the  next  lower  ply  to  make  a  single  matrix  stratum  of 
the  ccanbined  thicknesses.  A  laminate  of  N  plies  is  broken  into 
Ntl  matrix  strata  and  N  fiber  strata. 
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strata  Stiffness  Properties 


Each  matrix  stratum  is  a  homogeneous  layer  of  isotropic 
material.  Only  three  engineering  constants  are  required  to 
determine  the  stiffness  matrix  coefficients  of  an  isotropic 
homogeneous  plate  [Whitney:  12-13],  The  stiffness  matrix,  of 
matrix  strata  t,  [KmJt/  will  have  the  form: 


kj 
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1  Oy  1 

1  1^2 
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0  I 

1  1 
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1  ^  1 
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k3 
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I  1 

Jt 

1  0 
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0 
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(2) 


The  tensile  modulous  E,  Poisson's  ratio  v,  and  the  shear  modulous  G 
are  required  to  define  the  stiffness  matrix  coefficients  of  the 
matrix  strata.  G  is  customarily  determined  from  the  tensile 
modulous  and  Poisson's  ratio  by  the  relationship  G=E/2(l+v).  This 
reduces  the  required  number  of  constants  to  two,  E  and  v. 


Each  fiber  stratum  is  a  combination  of  matrix  and  fiber, 
dominated  by  the  fiber.  The  original  ply  thickness  (tp)  and  the 
fiber  radius  (r)  are  required  to  define  the  global  volume  fraction 
(Vf)  of  the  ply. 
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Figure  5.  Fiber  Strata  Pareuneters 


The  idealized  fiber  height  (hf)  of  each  stratum  is  calculated  by 
setting  the  fiber  volume  (nr^)  equal  to  the  idealized  fiber  height 
squared.  The  global  volume  fraction  (GVf)is  unchanged  and  the 
local,  fiber  stratum,  volume  fraction  (Vf)  can  be  defined  from  the 
fiber  radius  (r)  and  the  matrix  thickness  (t^). 


GVf 


Local  Vf 


xr 


(3) 


Each  fiber  stratum  is  modeled  as  a  homogeneous,  specially 
orthotropic  plate.  An  orthotropic  plate  requires  9  independent 
engineering  properties  En,  E22,  E33,  v^,  V13,  V23r  Gur  G13,  and  G23 
to  define  the  stiffness  matrix  coefficients.  A  representative 
volume  element  shows  that  Eii=E33,  vi3=vi2,  V23=V2i=(E22/Eii)*vi2,  «md 
Gi2=G23  which  reduces  the  number  of  required  constants  to  5. 
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fiber 

I  I  matrix 

1 

/ 

I 

612 

Figure  6.  Fiber  strata  Representative  Volume  Element 

The  modulus  and  Poisson's  ratio  in  each  direction  can  be 
determined  by  applying  the  fiber  and  matrix  moduli  and  Poisson's 
ratio  in  series  or  parallel  with  the  volume  fraction,  as  the 
geometry  dictates.  The  engineering  constants  are  determined  as 
follows: 


Ell 

-  VfEf  +  V„E„ 

Gi3 

=  VfGf  +  V„G„ 

t)l2 

- 

®22 

EfEm 

Gi2 

V23 

E22 

-  V12  ^ 

VfE„  +  V„Ef 

VfG„  +  V„Gi 

E33 

“  Ell 

G23 

-  Gi2 

^13 

-  VfUf  +  V„U„ 

Using  the  relationship  GsE/2(14v),  the  3-D  stiffness  matrix 

of  the  fiber  strata  in  each  ply  can  be  generated  if  the  constants 
E  and  V  are  known  for  the  matrix  and  the  fibers  involved. 
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The  on  fiber  axis  (x,y,z),  stiffness  matrix,  of  fiber  stratum 
r,  [Sf^,  will  have  the  specially  orthotropic  form: 


fSl 

93 

S3 

0 

0 

OJ 

f  1 

1  02  1 

IS3 

S2 

S4 

0 

0 

0  1 

I  £2  I 

103  1 

|S3 

S4 

Si 

0 

0 

0  1 

I  £3 1 

•i  X  ^ 

1  ^23  , 

^  0 

1 

0 

0 

S5 

0 

0  ^ 

1 

|Y23  , 

1  1 

1  1 

0 

0 

0 

0 

Sg 

o| 

Iv.j 

Jr 

[0 

0 

0 

0 

0 

IY12  J, 

(5) 


In  cases  when  the  fibers  are  not  aligned  with  the  primary  axis 
system  the  off-au^is  stiffnesses  will  be  required  to  define  the  ply 
stiffness.  The  off  fiber  axis  stiffness  [S],  may  be  calculated  by 
the  tensor  transformation  of  [S]'  as  follows: 

[S]  -  {Tr'[S]'{T}  .gv 


Using  the  transformation  matrix: 

n^  0  0  0  2  mn  1 

m^  0  0  0  -2  mn  I 
0  10  0  o' 

0  0  m  -n  0 

0  0  -t-n  m  0  I 

mn  0  0  0  m^  -  n^  J 

m  -  cos(6),  n  -  sin(6) 


{T} 


r  m" 

n^ 
0 
0 
0 

[-mn 


(7) 
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The  resultant  stiffness  matrix,  of  the  fiber  stratum  r,  will  be  of 


the  form: 


f  1 

f  Si 

S3 

SlO 

0 

0 

S5 

1 

f  "x  1 

1  Oy  1 

1  S3 

S2 

CO 

0 

0 

Se 

1 

I  By  I 

1  o,  1 

CD 

M 

0 

Sii 

Sl2 

0 

0 

Sl3 

1 

y 

1 

1 

I  e,  I 

^  X  ^ 

1  ”  1 

0 

0 

S7 

Sg 

0 

ly  I 

I  Yxx  I 

1  0 

0 

0 

Sg 

Sa 

0 

Jr 

1  S5 

Se 

^13 

0 

0 

S4 

1 

Jr 

lYjqr  Jr 

The  stiffness  matrices  of  both  the  fiber  and  matrix  strata 
Cein  be  fully  defined  if  the  fiber  radius  (r),  the  thickness  for 
each  ply  (tp),  the  global  volume  fraction  (Vf),  and  the  tensile 
modulus  (E)  and  Poisson's  ratio  (v)  of  the  fiber  and  matrix,  are 
known. 


Oeaeratioa  of  tbo  Ooveraing  Equations 

Using  Hamilton's  energy  principle,  the  variation  of  the 
volume  integral  of  the  internal  strain  energy  (II)  minus  the 
kinetic  energy  (T)  is  equal  to  the  work  of  the  applied  forces  (W). 

a  b  t/2  a  b 

//  /  (  6n  -  6T)dzdydx  -  J* j* 6Wdydx  (9) 

0  0  -t/2  0  0 
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For  simplicity,  this  document  only  examines  static  cases 
where  the  kinetic  energy  is  equal  to  zero.  The  internal  strain 
energy  is  broken  into  fiber  strata  and  matrix  strata 
contributions : 


n  *  +  Ilf 


(10) 


Internal  strain  energy  of  the  fiber  stratum  r  and  matrix 
stratum  t  are  functions  of  the  strain  vectors  and  stiffness 
matrices : 


(11) 


The  strain  vector  for  each  stratum  is  a  function  of  the 
displacements  u,  v,  and  w.  Only  the  linear  strain  terms  are  used 
and  the  thin  plate  assumption,  where  the  displacement  in  the  z 
direction  does  not  vary  in  z,  is  used: 


f  '^r,x  1 

\  1 

I  ''r,y  I 

1  Vt,y  1 

I  0  I 

1  0  1 

I  Ur,y  +  V,,,,  I 

,  +  W,,,  I 

|"t,*  +  I 

[  +  ^r,y  J 

(Vt,,  +  J 

The  total  strain  energy  is  the  summation  of  the  strain  energy 
of  all  the  strata.  M  is  the  number  of  matrix  strata  (N+l)  and  F 
is  the  number  of  fiber  strata  (N). 

Hi  -  2  "r  n.  -  2  n, 

r-1  t-1 

(13) 
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Work  is  a  function  of  the  load  and  displacement.  The  Strata 
Theory  assumes  small  displacements  and  only  considers  the 
displacements  in  the  z  direction  (w)  for  the  calculation  of  \#ork. 

W  =  p(x,y)*w(x,y) 

The  entire  expression  becomes  a  function  of  the  displacement 
functions  u,  v,  and  w  and  the  load  function  p(x,y).  The  first 
variation  can  be  taken  with  respect  to  the  displacements  (u,  v, 
and  w)  and  set  equal  to  zero  in  order  to  obtain  the  minimum  energy 
level  or  the  equilibrium  energy  state.  This  leaves  an  expression 
with  the  variation  of  the  derivatives  in  x,  y,  and  z  of  the 
displacement  functions.  Integration  by  parts  is  used  to  eliminate 
the  derivatives  of  variational  terms.  The  following  expression 
results  plus  a  set  of  boundary  conditions,  generated  from  the 
integration  by  parts. 
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1 


+  ^2t^t,xy  + 

k3t(Ut,yy  +  Ut,„  +  Vt,^  +  Wt,«  )}6Ut  + 

a  b  M  ^t  + 

1 

ffl*  f 

0  0  t-1  Zt- 

•k3t(Vt,*c  +  +  Wt,y,)}6Vt  + 

{Jt3t(Ut,«  +  +  Wt,yy  +  Wt,«)}SWt 

idzdydx  + 


[■{SirUr^xx  ®4r'^r,yy  ®6r^r,yy 

l(®3r  ■*■  ®4r)''^r,3y  "*■  ®5r(2^r,jiy  ^x,yx.)  "*■ 

l®7r(^r,22  ^r,X2  )  "*■  ^9r('^r,22  +  ''^r,y2  )}^^r  "*■ 


a  b  F  2r+ 


;/2*/ 


{®4r'^r,xx  ■*■  ®2r^r,yy  "*■  ®5r^r,xx 

(®3r  ■*■  ®4r)^r,:^y  "*■  ^erC'^rjyy  "*■  ^^r,xy  ) 

S8r(^r,22  +  Wr,y*)  +  Sg^(Ur,„  +  )}6Vr  + 


1 


idzdydx 


|{S7r(Wj.,«  +  Ur,«  )  +  | 

[S8r(''r,yy  +  '^r,yz)  +  +  “1,72  )}^x  \ 


a  b 

=  ff[p(x,y)*  dw]dydx 

0  0 


(15) 


The  above  equation  can  be  broken  into  3  equations  each  containing 
the  factor  of  6u,  6v  or  6w.  The  relationships  must  hold  for  an 

arbitrary  variation  of  unity,  so  the  variational  terms  may  be 
factored  out  leaving  the  following  set  of  differential  equations 
and  boundary  condition  equations  from  the  integration  by  parts. 
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r'- 


1/ 

a  b  f  M  I 

//  2  2i, 

Oor=lt-l|  y 

J 


^^3t('^t,yy  +  ^t,xy  "*■  ^t»2z  "*■  ''^t,X2  )  "*■ 
[k2tVt,xy  +  ^It^t.xx 

Slr'^r,xx  +  S4rUj.,yy  +  (83^  +  S^r)^r,xy 
S5r(2Ur  ,xy  +  Vr,xx  )  +  ®6r^r,yy  + 


Idz 


I  /^S5r(2Ur,xy  +  Vr,xx  )  +  ^6r^r,yy  + 

[^'-[s7r(Ur  ,22  ■*■  '’^r,x2  )  +  S9r(Vr  ,22  ^r,yz  )  J 


|dydx 


(16) 


a  b  F  M 

//  2  21, 

Oor*lt-lj  y 

I  ^ 

L 


f  ^3t(  ^t,xy  ■*■  '^tjXX  ■*■  Vt,22  ■*■  '*^t,y2  )  ] 

!/|k2tUt,xy  +  l^ltVt.yy 


Idz 


S2rVr,yy  +  (  +  S4r  )^r,xy  +  S4j.Vj.,xx  + 

S6r(2Vj-,xy  +  ^r,yy)  +  ^sj-Uj-^xx  + 

[S8r(  Vr  ,22  ^r,y2  )  +  S9r(Ur  ,22  ■*■  ^r,X2  )  J 


[dz 


dydx 


(17) 


a  b  F  M 

//2  2 

0  0  r«lt»l 


/ 


^3t('^t,x2  +  ^t,y2 
'*^t,xx  '*^t,yy  ) 


iz  + 


1 


I 


S7r(  '^r,xx  +  Ur,x2  )  + 
/•iS8r(  Vr,yz  +  ^r,yy)  +  l^dz 

■[s9r(Ur,y2  +  J 


a  b 

/  /  p(x,y  )lydx 


0  0 


J 


(18) 


Boundary  Conditions: 
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t-lr-l 


M  T 

‘-I  2 

t-lr-; 


r  r  *t+ 

/  (kitUt,x  +  lt2tVt.y}d2  + 


*t- 
*r  + 


/{sirUr,x  +  Sar^r.y  +  +  Vr,x)}dz 


dy  + 


r  2t+ 

/  {l^3tUt,y  +  k3tVt,x)dZ  + 


.«  2t- 


lo  I  Z 


L- 


/  {s5rUr,x  +  S6rVr,y  +  ^Ari^r,r  +  Vr,x)}dZ 

r- 

r[k3tut,z  +  ^3t^t,x 


[[^r<V 


*  +''r,x  )+Sr  (''r,y 


l*r+ 

^k.  j 


dx  + 

I 

J. 

1 

I  dxdy 


[ 


(19) 


r  1^  *t+ 

Hi 


r  2t+ 

•  +  J^3tVt:,x]<iz  + 


*t- 

*r  + 


I  I 

/{s5rUr,x  +  +  S4r(Ur,y  +  Vr,x)}dz  I 

L  *r-  J, 

r  *t+  1 

I  J*  {kjtUt,*  +  kitVt,x)d2  +  I 


J*{®3r'^r,x  ■*■  ®2r^r,y  +  S6r('^r,y  +  ^r,x 

r[k3tVt,,  +  k3tWt,y| 


dy  + 


-«  I  *t- 

J®  I  •r  + 


dx  + 


1-- 


l*t+ 


jf/oL  T., 

|^[®8r(^r,z  ■•■  ^r,y  )  +  S5j.(  Wj.^j5  +'^r,z)J^  j 


(20) 
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N  F 

t-lr-1 


r 

/  {k3tUt,z  +  k3tWt,xjdZ  + 

I  *"- 


I*  1 


dy  + 


1  jr{s9r(Vr,y  +  )}d2  . 

L  J 


J 

1 


r  2t+ 

/  {k3tVt,z  +  k3tWt,y}d2  + 

f  '  '  dx  + 

I  Zr. 

/{s8r(Vr,2  +  ''r,y)+  S9r(  ^r.z  +  Wr,x)}d2 


[jf  /r  [  ^  k3t«t,x 


i 


J 


(21) 


If  displacement  functions  can  be  found  that  satisfy  the 
differential  equations  and  the  boundary  conditions,  then 
compatibility  and  equilibrium  will  be  satisfied  and  the  strain  and 
stress  in  each  stratum  can  be  calculated. 


Strata  Model  Displacement  Assumptions 

The  assun^tions  of  the  strata  model  are  as  follows: 

1)  Only  linear  strain  terms  are  considered  significant. 

2)  The  displacement  in  w  does  not  vary  in  z  thus  ez=0. 

3)  There  is  constant  strain  through  the  thickness  of 
the  fiber  strata. 

The  main  premise  of  the  strata  model  is  that  the  fibers  and  thus 
the  fiber  strata  are  much  stiff er  than  the  matrix.  The  strata 
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model  assumes  constant  strain  through  the  thickness  of  the  fiber 
strata.  The  displacement  functions  in  x  and  y,  of  the  fiber 
strata,  will  not  vary  in  z.  Strain  in  the  matrix  strata  will  vary 
in  z.  The  displacement,  in  x  and  y,  of  each  position  in  the 
matrix  strata  will  be  a  function  of  its  z  coordinate. 

The  above  assumptions  will  generate  different  displacement 
functions  (u  and  v) ,  for  each  stratum.  To  maintain  continuity,  a 
continuous  displacement  function  must  be  maintained  through  the 
entire  thickness  of  the  plate. 

The  fiber  and  matrix  strata  interface  must  coincide  after 
deformation.  The  z  coordinate  of  the  bottom  of  the  matrix  strata 
above  each  fiber  strata  (Znibottcm) ,  and  the  top  of  the  matrix  strata 
below  the  fiber  strata  (2m(i+i)top)  must  be  mapped  to  the  fiber 
stratum  z  coordinate  Zfi. 


Undeformed  Deformed 


The  mapping  is  accon^lished  with  the  mapping  parameters  a  and  p. 
Each  matrix  strata  will  have  it's  own  a  and  p  mapping  parcuneters, 
that  are  a  function  of  the  surrounding  fiber  strata  coordinates 
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Zfi-1  and  Zf±.  The  linear  mapping  is  accon^lished  by  the  following 
equations . 


P%^itop  "  ^i-1 
Pinj_  ^m^bottcBi  ” 


(22) 


By  solving  both  equations  simultaneously  the  a  and  p  of  each 
matrix  stratum  can  be  calculated  as  follows. 


Pi 


Zf  -  Zf 

^1-1  M 
Zm^tOp  ~  Zgi^ 


f  ^1-1 

^1 

1  f  ^itop  * 

^ibottcn  1 

[  Zn^top 

^j^bottcra 

J  [  Zm^botten 

“  ^iBitcp  J 

(23) 


The  upper  most  and  lower  most  matrix  surfaces  will  be  mapped 
to  the  z  coordinate  of  the  free  surface  +/-  t/2.  In  the  above 
algorithm  this  gives  a  Zfo  of  t/2  and  a  Zf(r-n)  of  -t/2.  In  lay-ups 
with  all  plies  having  the  same  thickness  (tp),  and  idealized  fiber 
strata  height  (hf),  the  p  pareuneter  will  be  constant.  Matrix 
strata  with  positive  z  coordinates  will  have  negative  a  parameters 
and  matrix  strata  with  negative  z  coordinates  will  have  positive  a 

parameters.  The  linear  mapping  results  in  the  z  coordinate  for 
matrix  strata  t,  being  represented  with  (at+§t*z)  and  the  z 

coordinate  of  fiber  strata  r,  being  represented  with  the  constant 

To  generate  solutions,  the  Strata  Theory  assumes  that  the 
through-the-thickness  displacement  (in  the  z  direction),  is  a 
polyncmiial  function  of  z,  and  solves  for  the  displacement  in  x  and 
y.  The  constant  displacement  through  the  fiber  strata  can  be 
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modeled  with  a  polynomial  function  of  the  stratum  mid  plane 
distance  (Zf)  and  an  unknown  function  of  x  and  y.  The  displacement 

through  the  matrix  strata  varies  in  z  and  is  assumed  to  be  a 
polynomial  function  of  (a+pz)  and  the  same  unknown  function  of  x 

and  y.  The  displacement  in  the  z  direction  (w)  is  assumed  to  be  a 
function  of  x  and  y  only.  The  assumed  displacement  relationships 
are  represented  mathematically,  as  follows: 

N 

Uf(x,y,z)  =  y  zj  •  ii>i(x,y) 
iTo 

N 

Vf(x,y,z)  =  2  ^  ■  5i(x,y) 

i-O 

W£(x,y,z)  =  w(x,y) 

N 

Um(X/y/Z)  =  2(a«  +  •  ^i(x,y) 

iTo 

N 

v„(x,y,z)  -  •  5i(x,y) 

i-O 

WmCx^y^z)  =  w(x,y)  ^24) 

The  order  of  the  polynomial  (N)  can  vary  to  model  the 
complexity  of  the  loading  and  geometry.  More  complex  loading  and 
geometry  require  more  con^lex  displacement  functions.  The 
displacement  function  must  have  sufficient  order  to  suirvive  two 
derivatives  and  generate  complex  strain  functions.  A  polynomial 
order  of  5  has  been  found  to  have  sufficient  coii^>lexity  to 
represent  the  displacements  and  strains  of  the  simple  cases 
examined  in  this  document  and  shall  be  used  in  all  the  example 
cases  unless  otherwise  stated. 
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Solution  Procedures 


Using  the  above  displacement  functions  to  derive  the  strain 
energy  vectors,  the  2  differentials  in  the  governing  equations  can 
be  evaluated  before  the  first  variation.  The  governing  equations 
become  functions  of  z,  the  functions  ^(x,y),  w(x,y)  and  the 

load  function  p(x,y).  The  first  variation  is  taken  on  the  W,  ^  and 
w  functions,  generating  the  Strata  Theory  governing  equations  as 
follows : 


a  b  N  F  M 

fflll 

0  0  i-Or-lt-l 


(k2t  +  ^3t  )^i,xy^£  +  | 

(S3r  +  S4r)Ci,xyZ£  +  S4r^i,yyZ£  +  jdyd*  =  0 


Slr^i,xxZ£  +  rXxZf  + 

2S5r^i,xx^f  ®6r^i,xx^f 


J 


(25) 


a  b  N  F  M 

//  2  S  2 

0  0  i-Or-lt-1 


^3t  ( ^iZin,z  ^,yZm,2  )  "*■  ^ 

(^2t  ■*"  ^3t  )^i,xy^f  ^lt^i,yyZin 
(®3r  ■*"  ®4r  )^i,xyZf  +  ®4r^i,xxZf  +  Idydx  «  0 
S2r^i,yyZ£  +  S5rH»i,3cxZf  + 

[®6r(  ■'Pi.yyZf  +  2^i^xyZf  ) 


(26) 


23 


1 


r^3t('"^,xx  ^,yy )  "*■  1 

a  b  N  F  M  I  I  a  b 

//  E  2  2  ^3t(^i,x  +  Ci,y)Zro,z  +  dydx  -  JJ p(x,y)ciydx 
0  0  i-Or»lt-ll  I  0  0 

[®7r''^,xx  +  S3j-W^yy  +  2S9j.VryyJ 


where 


*t+ 

Zm  “  /(«t  +  M)^dz 

*t- 


J*  Z^rdZ 


*r- 


(27) 


(28) 


and 


Zm  -  /(at  +  PtZlf'"'"^’dz 

Zt- 

Zn>,2  -  /i  *  (at  +  Pt2/"’"^^dz 


Zi  -  7z<r^>dz 


"■frZ 


/i  *  4^  ^>dz 


Zm,z  -  /-ij  *  Pt(at  +  z|,2  -  / -ij  *  4r^^'^^dz 

2t-  z*. 

j  =  0..  N  ^29) 

The  derivation  also  generates  the  following  set  of  boundairy 
conditions  from  the  integration  by  parts: 


[ 


M  F  N 

111 

t-lr-li-0 


r  2  2  2  * 

I  J^lt^i,xZ«+J^2t^i,yZm  +1 

0  6tpi  dy  + 

[®3r^i,yZ£  +  S5r(  +  ^i,x  )Zf  J 


1 


So  ^i 


J^3t(^l,y  +  ^3t^i,x  )Zm  +  Ssr^ljxZf  +  "" 


S6r^i,yZf  +  S4r(  ^i,y  +  ^i,x  )Zf 


dx 


Jo 


(30) 
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,  b  r^3t^i,y2^  +  +  Ssr^ijX^f  +1  ' 

1 1  8^ •  dy  +  1 

«  r  N  1*^®  S4r(ll>i,y  +  ^i,^)Zi  +  85^Ci,yZi  I 

2  2  2  I  b  I 

t.lr-li-0|  +  S3,ti>i,,z2  +  , 

[•’”  ^  t  S.rC'Pl,,  *  )Z|  +  S,,Ci,,2|  J^  J 


(31) 


N  F  N  [  k, .  11) ,  Z  +  k W  Z  +  S,  W  Z,  +  S„  W  Z.  T( 

"  ^  "  I  Jo  |_  3t^i  m,t  3t  ,x  m  7r  ,x  f  9r  ,x  f  Jj 

2  2  Z  I  ^  r  -fc 

t-lr-li-O  r  ^  _  +k  W  Z  +S-  W  Z  +  S_  W  Z,  fd 

Ijo  L  3t^i  m,*  3t  ,y  m  9r  ,x  f  8r  ,y  £  Jj 


dy  + 


(32) 


With  the  correct  choice  of  the  ij)(x,y),  ^(x,y)  and  w(x,y) 
functions  and  the  load  functions  p(x,y),  the  x  and  y  dependence 
can  be  evaluated  directly  and  factored  or  integrated  out  of  the 
governing  equations.  The  governing  equations  can  then  be  solved 
and  the  displacement  functions  determined. 
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III.  strata  Model  For  a  Simply  Supported  H<mogeaeou8  Plate 


Simply  Supported  Homogeneous  Plate  Equations 


For  verification  of  the  theory,  the  strata  model  can  be 
solved  for  the  simple  case  of  a  singly  supported  homogeneous 
plate.  For  a  homogeneous  plate  the  fiber  strata  stiffness 
constants  (Sf)  reduce  to  the  form  of  the  matrix  strata  coefficients 
(]^)  where  Si=S2=  ki,  S3=  ka,  S4=S7=S8=  ka  and  S5=S6=S9=  0.  The 
governing  equations  25-27  reduce  to: 


a  b 

M  N 

r  kst  (  2  Ipi^yyZj,,  + 

!  +  )  + 

1 

// 

I  2 

o 

.  ^  y 

jdydx  >  0 

0  0 

t- li  -C 

•[(*^2t  +  k3t  )2Cj^^ 

J 

a  b 

M  N 

rJ^3t(2^i,xxZm  + 

1 

// 

2  1 

1  , 

O 

,  ^  c-  y 

jdydx  -  0 

0  0 

li.  *( 

>[(  ^2t  +  ^3t  )2^1 

.,xy^in 

*■  1  2  ,yy2in 

J 

(33) 


(34) 


ab  M  N  fk3t(2w^3nj  +  2w^yy)  +1 


a  b 


t^li  ■o[  ^3t  (  ^i,x  ^i,y  )^in,z  J 


jdydx  ~  ff  p(x,y)dydx 
0  0 


(35) 


Equation  28  must  be  modified  to  allow  displacement  to  vary  in  z 
through  the  entire  thickness.  Equation  28  becomes: 
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*t+ 

Z„  -  J* 

*t- 


The  boundary  conditions  reduce  to: 


(36) 


^  r  I  ^ 

'■  °[JC  SViLkatSCifi.,  +  katCi,.  )Zi  J 

•■  "  [jf6ti[*'3t2('l>i,,zi  +  k3rfl.xZi)T‘’y  ■"! 

2  2  b  “  ° 

t.u.olj^8^l[^^t2Vi,X  *  I'ltZSiXl*'  j 

M  »  [jf  »«[k3.  (  fiZ...  *  2«,A 

’I"**  j 


(37) 


(38) 


(39) 


Functions  to  represent  the  x  and  y  dependence  of  the 
displacement  are  required  to  solve  the  homogeneous  plate  governing 
equations  (33-39).  Using  proper  Sine  and  Cosine  series  functions 
for  the  X.  and  load  functions,  can  force  the  x  and  y  dependence 

to  become  constant  throughout  each  of  the  governing  equations. 

The  displacement  and  load  relationships  [2:  30]: 
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i|>i(x,y)  -  Ai  •  Cos(px)Sin(qy) 
^i(X/y)  -  Bi  •  Sin(px)Cos(qy) 
''(X/Y)  “  C  •  Sin(  px  )Sin(  qy  ) 
p(x,y)  -  P  •  Sin(px)Sin(qy ) 

i  -  0, 1,. . N 


(40) 


generate  x  and  y  dependent  derivatives  in  equations  33-35  that  go 
to  zero  over  the  integral  of  the  plate.  The  x  and  y  dependence  is 
constant  in  all  the  terms  of  each  governing  equation  and  can  be 
factored  out,  leaving  a  set  of  governing  equations  as  a  function 
of  z  and  the  coefficients  Ai,  Bi  and  C.  The  z  integrals  can  be 
evaluated  directly,  leaving  the  governing  equations  as  function  of 
the  coefficients  Ai,  Bi  and  C. 

The  load  p(x,y)  in  equation  40,  represents  a  distributed 
pressure  that  varies  in  x  and  y  with  a  maximum  pressure  at  the 
center  of  the  plate.  For  a  plate  singly  supported  on  all  four 
edges,  this  should  create  a  domed  displacement  function  with  the 
mciximum  deflection  at  the  center.  More  complex  loading  can  be 
modeled  with  a  full  Fourier  Sine  series  [6:  549]. 

The  governing  equations  can  be  transformed  into  a  linear 
system  of  equations  as  follows: 


(Ai 

[M]^Bi 

[  c 


i  -  0..N 


(41) 


M  will  be  [2*(N+1)+1]  X  [2*(N+1)+1]  matrix  while  the  other  two 
matrices  will  be  2*(N+1)  x  1  vectors.  Guassian  elimination  can  be 
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used  to  solve  for  the  coefficients  Ai,  Bi,  and  C.  The  coefficients 
are  then  used  to  generate  the  displacements,  strains  and  stress 
throughout  the  laminate. 

The  boundary  conditions  for  a  singly  supported  plate  have  v, 
w  and  Ox  equal  zero  a'  x=0  and  x=a,  and  u,  w,  and  Oy  equal  to  zero 

at  y=0  and  y=b.  This  satisfies  boundary  conditionof  equation  39. 
By  definition  6u=0  at  y=0,  and  y=b  and  6v=0  at  x=0  and  x=a,  which 

satisfies  boundary  conditions  of  equations  37  and  38. 


Simply  Supported  Homogeneous  Plate  Example 


A  homogeneous  plate  test  case  was  solved  using  the  above 
procedure.  The  following  pareuneter  are  required  to  define  the 
plate. 


Table  1.  Homogeneous  Plate  Parameters 

Parameter _ 

Number  of  plies  (n) 

Ply  thickness  (h) 

Laminate  thickness  (t) 

Plate  length  (a) 

Plate  width  (b) 

Tensile  Modulus (E) 

Poisson's  ratio  (v) 

Shear  Modulus  G=E/2(l+v) 

Polynomial  Order  (N) 


Value 

2 

1/2  in 
1  in 
10  in 
10  in 

3.0*10®  psi 
0.4 

1.071*10®  psi 
5 


With  the  above  parameters  the  stiffness  and  displacements  of 
a  homogeneous  plate  can  be  fully  defined  and  the  system  of 
equations  solved  for  the  x  and  y  dependence  coefficients.  The 
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laminate  thickness  (t)  and  the  load  (P)  were  chosen  to  normalize 
the  plate  dimensions  and  displacements  to  the  thickness  and  to 
normalize  the  plate  moduli  to  the  load.  The  stratified  plate 
model  displacements,  strains  and  stresses,  of  a  singly  supported 
homogeneous  plate  are  shown  in  Figures  8  through  22.  The  x  and  y 
dependence  has  been  factored  out  to  give  the  through  the  thickness 
distributions.  All  the  length  dimensions  have  been  normalized  to 
the  plate  thickness  (h)  and  the  modulous  dimensions  have  been 
normalized  to  the  load  (P). 

3.4 


(X10--6) 

-3.4 

-0.5  0  0.5 

Normalized  z 

Figure  8.  U(z)  Displacement  for  a  Homogeneous  Plate 
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-0.5 


0 


.5 


Normalizsd  z 

Figure  10.  W(z)  Displacement  for  a  Homogeneous  Plate 
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-0.5 


0 


0.5 


Normalized  z 

Figure  12.  Epsilon  y  (sy)  Strain  for  a  Homogeneous  Plate 
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Normalized  z 


Figure  17.  Sigma  x  (ax)  for  a  Homogeneous  Plate 


Normalized  z 


Figure  18.  Sigma  y  (c^)  for  a  Homogeneous  Plate 


Normalized  z 

Figure  19.  Sigma  z  (oz)  for  a  Homogeneous  Plate 


Normalized  z 


Figure  20.  Tau  xy  (xxy)  for  a  Homogeneous  Plate 


Normalized  z 


Figure  21.  Tau  xz  (xxz)  for  a  Homogeneous  Plate 


Normalized  z 

Figure  22.  Tau  yz  (xyz)  for  a  Homogeneous  Plate 
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simply  Supported  Homogeneous  Plate  Discussion 

The  results  shown  in  figures  8-13  are  the  expected  linear 
displacement  and  normal  strain  functions  for  a  homogeneous  plate. 
The  upper  half  of  the  plate  is  in  tension  from  the  bending,  but 
experiences  a  negative  displacement  due  to  the  rotation  allowed  by 
the  B.C..  The  constant  displacement  in  z  is  a  result  of  the 
strata  model  displacement  assumptions  and  corresponds  with  most 
linear,  small  displacement  plate  models.  The  parabolic  stress  x^z 
and  Tyz  strain  and  stress  distribution  is  a  result  of  mid  plane 

carrying  the  maximum  transverse  shear,  and  the  transverse  shear 
going  to  zero  on  the  top  and  bottom  surfaces.  The  existence  of  a* 
is  a  result  of  the  e*,  £y  and  terms  of  the  strain  vector 
multiplying  out  with  terms  in  the  stiff ess  matrix.  The  stress 
distribution  can  only  be  explained  by  assuming  that  half  the  load 
is  applied  to  the  lower  surface  and  half  the  load  is  applied  to 
the  upper  surface. 
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IV. 


Solution  to  a  Simply  Supported  Plate  with  Cross  Plies 


Composite  Plate  with  Cross  Plies  Equations 

For  a  siii5)ly  supported  square  conqposite  plate  with  only 
90  degree  plies,  the  fiber  strata  stiffness  constants  S5,  Ss, 
are  equal  to  zero  and  equations  25-27  reduce  to: 


r  k3 ^ ^ i ^yyZjn  +  ^ 

//  i  i  i  +  k3t)ti,xyZ?  +  +  idyd* 

0  0  r=lt*li*0  2  2  2 

[(®3r  +  ®4r)^i,xy^f  ®4r^i,yy^f  Slr'’l’x.xx2f 


ab  F  M  N  I  I 

SSI  II\^'^2t  ^3t  )^i,xy^f  jdydx 

Oor-lt-li-0  2  V  2  v  21 

[(^3r  +  ®4r  )^i,xy^f  S4rCi,xx2f  ®2r^i,yy^f  J 


r^3t('^,xx  '^.yy  )  1 

A  b  F  M  N  s  Id 

ff  I  I  i!k3t(^i,x  +  ^i,y)Zm,z  +  -  //p(x,y)dydx 

0  0  r-lt-li-ol  '  0  0 

L^Tr'^/Xx  +  ®8r'*^,yy  J 


*t+ 


Zm  “  /(«t  + 


f 


»t- 


*r- 


0  or 
and  S9 


0 

(42) 

0 

(43) 


(44) 

(45) 
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The  boundary  condition  equations  20-32  reduce  to: 


222] 

t-lr-li-O 


/o  ^i 


r,  2  .  ^  2  1 

r  1 

1  ^lt^i,x^m  ■*'^2t^i,y^m  "*■ 

|^®lr^i,x^f  ■*‘®3r^i,y^f 

dy  + 

0  1 

■  1^3t(3l>i,y  +  k3t^i,x)Z^  +  ■ 

b  1 

S4r(ll)i^y  +  )Z| 

dx  1 

0  J 

(46) 


r 


M  F  N 

22  2i 

t-lr-li-0| 


V  2  * 


/o^i 


®3r^i,x2i  ■*■  ®2j^i,yZf 


0 

b 

dx 


J 


(47) 


H  F 


»  +  ,Zj  Jdy  +1 

t-lr-li-O  |-  ^  +k  w  Z  +S-W  Z.  dx  I 

l^jo  L  3t^i  m,*  3t  ,y  m  8r  ,y  £  J 


(48) 


The  boundary  conditions  for  a  singly  supported  plate  have  v, 
w  and  Ox  equal  zero  at  x=0  and  x=a,  and  u,  w,  and  Oy  equal  to  zero 

at  y=0  and  y=b.  This  satisfies  boundary  condition  of  equation  48. 
By  the  simply  supported  definition  6u=0  at  y=0,  and  y=b  and  6v=0  at 

x*0  and  x=a,  which  satisfies  boundary  conditions  of  equations  46 
and  47. 
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The  Tj),  w  and  p  functions  of  equation  40  can  be  used  in 
equations  42-45.  The  x  and  y  dependence  derivatives  once  again  go 
to  zero  but  are  constant  in  each  equation  and  can  be  factored  out. 
The  z  integral  can  be  evaluated  directly,  leaving  the  equations  as 
a  function  of  Ai,  Bi  and  C.  The  system  of  equations,  can  then  be 
solved  for  the  Ai,  Bi  and  C  coefficients. 


Simply  Supported  Square  Plate  with  0/0  Plies  Example 


The  exan^le  plate  in  table  1  is  used  for  the  con^xjsitp  plate 
test  case.  Additional  composite  parameters  must  be  supplied  to 
fully  define  the  composite  plate. 

Table  2.  Con^site  Plate  Parameters 


Parameter 

Value 

Number  of  plies  (n) 

2 

Ply  thickness  (h) 

1/2  in 

Laminate  thickness 

1  in 

Plate  length  (a) 

10  in 

Plate  width  (b) 

10  in 

Ply  orientation  angles  6i 

0  radians 

62 

0  radians 

Matrix  Tensile  Nodulous (Em) 

3.0*10®  psi 

Matrix  Poisson's  ratio  (Vn) 

0.4 

Matrix  Shear  Nodulous  Gin=E/2  ( 1+Vn») 

1.071*10®  psi 

Fiber  Tensile  Modulous(E£) 

2.0*10'^  psi 

Fiber  Poisson's  ratio  (Vf) 

0.35 

Fiber  Shear  Nodulous  Gf=Ef/2(l+Vf) 

1.407*107  psi 

Global  Fiber  Volume  Fraction 

0.5 

Global  Matrix  Volume  Fraction 

0.5 

Fiber  radius 

•Jl  /  A-ja  in 
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The  ply  thickness  was  chosen  to  generate  a  total  thickness  of 
1,  to  automatically  normalize  the  results  to  the  thickness.  The 
fiber  radius  (not  realistic)  was  derived  from  a  global  volume 
fraction  of  1/2.  From  the  supplied  information  in  Table  2  the 
required  strata  model  parameters  may  be  calculated  as  follows: 


Strata  Model  Parameters 

Parameter 

Value 

Polynomial  Order  (N) 

5 

Number  of  Matrix  Strata  (n+1) 

3 

Number  of  Fiber  Strata  (n) 

2 

Idealized  fiber  height  (b) 

V2  /  4 

Matrix  thickness  (tm) 

0.  5  -  V2  /  4 

Local  Fiber  Strata  Volume  Fraction 

V2  /  2 

Local  Fiber  Strata  Volume  Matrix 

1  -  J2  /  2 

Mapping  Parameter  Pi 

2  +  Jt. 

Mapping  Parameter 

-(1  +  J2)  /  2 

02 

0 

03 

(1  +  V2)  /  2 

p  (x/a) 

ji/10 

q  (x/b) 

ji/10 

The  idealized  composites  plate  is  is  pictured  as  shown  below: 


X 


Figure  23.  Exan^le  Con^site  Plate. 
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displacements,  strains  and  stresses  of  a  simply  support 
plate  with  0/0  degree  plies  are  shown  in  figures  23-37 
y  dependence  has  been  factored  out. 


-0.5  0  0.5 

Normalized  z 

24.  U(z)  Displacement  for  a  0/0  Ply  Plate 


Normalized  z 
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Figure  25.  V(z)  Displacement  for  a  0/0  Ply  Plate 
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Figure  26.  W(z)  Displacement  for  a  0/0  Ply  Plate 


Normalized  z 


Figure  27.  Epsilon  x  (ex)  Strain  for  a  0/0  Ply  Plate 


Normalized  z 


Figure  28.  Epsilon  y  (ey)  Strain  for  a  0/0  Ply  Plate 
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Figure  29.  Epsilon  z  (8z)  Strain  for  a  0/0  Ply  Plate 
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Gaama  yz 
(xl0"-6) 


Normalized  z 


Figure  32.  Gamma  yz  (yyz)  for  a  0/0  Ply  Plate 


Normalized  z 


Figure  33.  Sigma  x  (ax)  for  a  0/0  Ply  Plate 


Normalized  z 


Figure  34.  Sigma  y  (oy)  for  a  0/0  Ply  Plate 


Sigaa 


Normalized  z 


Figure  35.  Sigma  z  (oz)  for  a  0/0  Ply  Plate 


8 


Normalized  z 


Figure  36.  Tau  xy  (xxy)  for  a  0/0  Plate 


Normalized  z 


Figure  37.  Tau  xz  (xxz)  for  a  0/0  Plate 
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-0.5  0  0.5 

Normalized  z 

Figure  38.  Tau  yz  (xyz)  for  a  0/0  Plate 

Simply  Supported  0/0  Composite  Plate  Discussion 

The  displaceinents  shovm  in  figures  21-23  are  expected  for  a 
con^osite  plate.  The  outer  matrix  strata  exhibit  greater 
displacements,  in  x  and  y,  than  the  stiff er  fiber  strata.  The 
ductile  matrix's  displacement  increases  with  the  distance  from  the 
mid  plane.  The  maximum  displacement  is  in  the  z  direction  as  can 
be  expected  from  the  loading.  The  m£ucimum  strain  is  found  to  be 
f^he  Yjv  shear  strain  from,  the  coupling  of  the  normal  strains.  The 
shear  strains,  and  Yyz  are  discontinuous  functions  because  the 
strata  model  assun^ion  that  the  fiber  strata  strain  does  not  vary 
in  z.  This  creates  a  plain  strain  assumption  in  z  of  the  fiber 
strata.  Because  of  this  anomaly  the  Yxi  and  Yys  strains  from  the 
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strata  model  should  not  be  relied  on.  The  existence  of  Oz  is  a 
result  of  the  ex,  ey  and  Yj^  terms  of  the  strain  vector  multiplying 
out  with  the  Sio,  Sa  and  s^  terms  of  the  fiber  strata  stiff  ess 
matrix  and  the  k2  and  ka  terms  of  the  matrix  strata  stiffness 
matrix.  The  stress  distribution  can  only  be  explained  by  assuming 
that  half  the  load  is  applied  to  the  lower  surface  and  half  the 
load  is  applied  to  the  upper  surface.  The  model  produces  fairly 
significant  normal  stresses  in  z.  The  transverse  normal  stress, 

Oz,  is  on  the  same  order  of  magnitude  as  the  shear  stress  This 

a  peculiarity  of  the  loading  assus^tions  in  the  strata  model.  The 
presence  of  the  stiffer  fibers  creates  a  discontinuous  stress 
distribution  in  all  directions.  The  maximum  stress  is  experienced 
by  the  fiber  strata  while  the  maximum  normal  strain  is  experienced 
by  the  outer  matrix  strata. 

The  stress  and  strain  distributions  generated  by  the  strata 
model  parallel  those  generated  by  a  Classic  Laminated  Plate  model 
analysis  of  a  square  plate  under  a  uniform  moment  about  the  x  and 
in  y  axes.  The  Ex,  Cy  and  strains  and  the  ox,  Oy  and  stresses 

generated  by  a  Classic  Leuninated  Plate  model  analysis  with  a  plane 
stress  assumption,  exhibit  the  same  patterns  of  stress  and  strain 
as  the  strata  model. 

Square  Plate  vritb  0/90  Degree  Plies  Example 

The  displacements,  strains  and  stresses  for  a  0/90  plate  are 
shown  in  figures  38-52.  The  x  and  y  dependence  has  been  factored 
out  to  give  the  through  the  thickness  distributions. 
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Figure  40.  V(z)  Displacement  for  a  0/90  Laminate  Plate 
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Figure  42.  Epsilon  x  (ex)  Strain  for  a  0/90  Laminate  Plate 
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Figure  53.  Tau  yz  {xyz)  for  a  0/90  Laminate  Plate 

Simply  Supported  0/90  Composite  Plate  Discussion 

As  expected,  a  larger  displacement  and  normal  strain  in  x  is 
experienced  by  the  90  degree  ply,  because  of  the  orientation  of 
the  fibers.  The  shear  strain  is  equal  in  both  plies  because  of 

the  counter  acting  effects  of  the  fiber  strata.  The  shear 
strains,  Yxz  and  Yyz  sre  discontinuous  functions  because  of  the 

strata  model  assumption  that  the  fiber  strata  strain  does  not  vary 
in  z.  The  maximum  strain  is  experienced  by  the  fiber  strata  in 
the  transverse  shear  direction  (Yj^)*  Discontinuous  stress 

distributions  are  experienced  in  both  plies  because  of  the 
presence  of  the  stiff er  fiber  strata.  The  maximum  stress  is 
experienced  by  the  fibers  in  the  transverse  shear  direction, 
because  of  the  plain  strain  condition. 


59 


V.  Cylindrical  Bending 


The  singly  supported  boundary  conditions  of  sections  III  and 
IV  are  only  useful  for  plates  supported  on  all  four  edges.  A 
plate  supported  on  two  opposing  edges  exhibits  cylindrical  bending 
where  displacements  and  strains  in  one  directions  go  to  zero. 

This  type  of  boundary  condition  can  be  modeled  with  the  strata 
theory  by  letting  the  y  dimension  (b)  go  to  infinity  and  examining 
the  plate  at  the  center  [2:  57-58].  Letting  b  go  to  infinity 
makes  the  displacement  v  and  derivatives  with  respect  to  y  go  to 
zero.  This  simplifies  the  governing  equations  and  produces  zero 
values  for  the  displacement  v,  and  the  strains  ey,t^  and 


Cylindrical  Bending  Governing  Equations 

The  governing  equations  (25-27)  for  a  plate  with  cross  plies 
in  cylindrical  bending  reduce  to: 


i  i.Or-lt-l^SirH>i,xxZf  +  +2S5j.ll>i^3„z|  j 


(49) 


a  N  F  M  a 

/  X  ^  ^3t^i,x2in,z  ^  ^7r^,xx  ]ix-/p(x)dx  (50) 

0  i*«0r«lt-l  0 
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Z, 


*t+ 

®  “  /(“t  + 

at- 


Zf 


ar  + 

f 


Zr- 


(51) 


The  boundary  Conditions  30-32  reduce  to 


H  F  N  rb 


2;  2 ^]S  ^ I dy J 


(52) 


M  r  N  r  r  Ta 

^  ^  k,.  1I)  .Z  +  k,^w  Z  +  s_  w  Z.  +  w  Z^  I 

"  "  I  3t^i  ]n,z  3t  ,x  m  7r  ,x  f  9r  ,x  f  L 

t-lr-li-QL*-  JO 


(53) 


The  tp,  w  and  p  functions  of  equation  40  can  be  used  if  the 
y  dependence  is  deleted.  The  x  and  y  dependence  takes  on  the 
following  form 


^i(Xf  y)  -  •  Cos(px) 

y )  •=  C  •  Sin(  px) 
p(x,y)  -  P  •  Sin(px) 


Using  the  above  relationship  in  equations  49-51,  allows  the  x  and 
y  dependence  derivatives  to  be  factored  out,  and  the  system  of 
equations  solved  for  the  Ai  and  C  coefficients. 

The  boundary  conditions  for  a  plate  in  cylindrical  bending 
have  V,  w  and  a*  equal  to  zero  at  x=0  and  x=a.  This  satisfies 

boundary  condition  53.  A  singly  supported  plate  has  zero 
curvature  at  x=0  and  x=a,  so  the  and  Ci,x  terms  in  equation  52 
go  to  zero  satisfying  the  boundary  condition. 
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0/0  C<MpoBite  Plate  in  Cylindrical  Bending  Example 

The  same  san^le  plate  described  in  Table  2  and  3  is  used  for 
cylindrical  bending  accept  that  the  y  dimension  (b)  is  set  at 
infinity.  The  displacements,  strains  and  stresses  of  a  con^site 
plate  with  0/0  plies  in  cylindrical  bending  are  presented  in 
figures  53-60.  The  x  and  y  dependence  has  been  factored  out  to 
give  the  through  the  thickness  distributions. 


Figure  54.  U(2)  for  a  0/0  Plate  in  Cylindrical  Bending 
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sigma  x 


Figure  58. 


Normalized  z 


ax(z)  for  a  0/0  Plate  in  Cylindrical  Bending 
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Tau 


-0.5 


0 


0.5 


Normalized  z 


Figure  61.  for  a  0/0  Plate  in  Cylindrical  Bending 


Cylindrical  Bendijig  0/0  Plate  Diacuaaion 

The  displacements  shown  in  figures  52-53  are  expected  for  a 
singly  supported  con^site  plate  in  cylindrical  bending.  The 
outer  matrix  strata  exhibit  greater  displacements  than  the  stiffer 
fiber  strata.  The  ductile  matrix’s,  displacement  increases  with 
the  distance  from  the  mid  plane.  The  presents  of  the  stiffer 
fibers  creates  a  discontinuous  stress  distribution  in  all 
directions.  The  plate  still  experiences  stress  in  the  y  direction 
from  a  Poisson's  contraction  effect.  The  maximum  normal  stress  is 
experienced  by  the  fiber  strata  vdiile  the  maximum  normal  strain  is 
experienced  by  the  outer  matrix  strata. 


66 


0/90  Composite  Plate  in  Cyliadrical  Bending  Example 

The  same  san^le  plate  described  in  Table  2  and  3  is  used  for 
cylindrical  bending  accept  that  the  y  dimension  (b)  is  set  at 
infinity  and  02  is  set  a  90  degrees.  The  displacements,  strains 

and  stresses,  independent  of  x  and  y,  of  a  coo^site  plate  with 
0/90  plies  in  cylindrical  bending  are  presented  in  figures  61-68. 


Normalized  z 

Figure  62.  U(z)  for  a  0/90  Plate  in  Cylindrical  Bending 


67 


68 


Normalized  z 


Figure  66.  o^iz)  for  a  0/90  Plate  in  Cylindrical  Bending 
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Normalized  z 


CT-f 


) 


for  a  0/90  Plate  in  Cylindrical  Bei 
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Tau  xz  0 


-  9  -I - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 

-0.5  0  0.5 

Normalized  z 

Figure  69.  Txz(.’  for  a  0/90  Plate  in  Cylindrical  Bending 

Cylindrical  Beading  0/90  Plate  Discussion 

As  expected,  a  larger  normal  displacement  and  strain  in 
X  is  experienced  by  the  90  degree  ply,  because  of  the  orientation 
of  the  fibers.  The  reverse  is  true  for  the  0  degree  ply  in  the  y 
direction.  Once  again  we  get  discontinuous  Ym  because  of  the  zero 

strain  gradient  assuiq>tion  in  z.  Discontinuous  stress 
distributions  are  experienced  in  both  plies  because  of  the 
presence  of  the  stiffer  fiber  strata.  The  maximum  stress  is 
experienced  by  the  fibers  in  the  transverse  shear  direction. 
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VI.  Rxtz  Polynomial  Solution  to  Composites  with  Angle 

Plies 

With  angle  plies  included  in  the  governing  equations  the 
fiber  strata  stiffness  matrix  coefficients  Sj,  Sg  and  Sg  remain  in 
the  governing  equations.  With  the  Sine/Cosine  V  and  ^  functions, 
the  X  and  y  dependent  derivatives  are  not  constant  in  each  of  the 
equations  and  can  not  be  factored  out.  Complete  functions  to 
represent  the  x  and  y  dependence  of  the  displacement  are  required 
to  solve  the  governing  equations.  A  Ritz  polynomial  series  can  be 
used  for  the  x  and  y  dependence  functions  W  and  The  x  and  y 

dependent  derivatives  can  then  be  evaluated  and  integrated  out  of 
the  governing  equations  and  the  equations  solved.  The  Ritz 
displacement  relationships  are  represented  as  follows 

N  H 

i-lj-l 
19  H 

i-lj-1 
N  H 

Wi,j(x»y)  -  2 

i-lj-i  (55) 

The  double  series  can  be  transformed  into  a  single  series  of  the 
form: 
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a 

^i(x,y)  -  2)  Ai  •  fi(x,y) 


( X ,  y ) 


Wi(x,y)  = 


P(x,y) 


i-l 

2)  Ci  •  fi(x,y  ) 

i-l 


(56) 


Different  applied  loads  p(x,y)  can  also  be  represented  by  a 
Ritz  polynomial  series.  A  single  constant  load  term  independent 
of  X  and  y  represents  a  constant  distributed  pressure. 

For  the  Ritz  polynomial  displacement  assumptions  the 
governing  equations  (25-28)  become: 


^3 1  ( ^ i ,yy2in  "*■  ^ 

|<>^2t  +  J^3t)^i,xyZf  +  +  ) 

j(S3r  +  S4r)^i,xyZf  +  S4r^i,yyZf  +  j^Y^X 

|Sir^i,xxZf  +  Ssr^i^xxZf  +  | 

[2S5rT)>i,xxZf  +  J 


a  b  N*  F  M 

fflll 

0  0  i=lr  =  lt-l 


a  b  F  M 

//2  2  2 

0  0 


|^^2t  ^3t  )^i,xyZf  ■*■  ^lt^i,yyZm  I 


|^^2t  ^3t  )^i,xyZf  ■*■  ^lt^i,yy2m 

j(S3r  +  S4r  )^i,xyZf  +  S4j.^i^xxZf  + 

_2  2 
|S2r^i,yyZ£  +  Ssr^i  ,xxZf  + 

®6r(^i,yyZf  +  2^i^jjyZf  ) 


Idydx 


(57) 


(58) 
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1 


J/iii 

0  0  i*lr=lt«l 

®7r''^i,xx  +  Sgj-Wi^yy  +  2S9j.Wi^yy  j 


a  b 

dydx  -  JJ  p(x,y)dydx 
0  0 


where 


5m  =  /(«t  +  M)^dz 


f  z^dz 


The  boundary  conditions  of  equations  30-32,  for  a  Ritz  Polynomial, 
become: 


M  r 

222 

t-lr-li-l 


[  r ,  2  ,  .  2  2  -1®  1 

I  >  I  +  J^2t^i,yZm  +  +  I  , 

Jo  ^i  I  I 

I  [®3r^i,yZi  +  S5r(^i,y  +  ^i,x  )2f 

!  a  ^3t(^i,y  ®5r^i,x^f  "*■ 

'  fn  ^  i  5  7  ^ 

[^er^iA  ^  «4r(^i,y  +  ^i,x  )Zi  Jo  J 


(61) 


M  F 

22  2 

t-lr-li-l 


\  J>  J^3t^i,yZm  +  ^  S5rtl,xZf  +  ' 

r°  "■  [S4r(^i,y  +  Ci,^)Zi  +  S6r^i,yZf  ]„ 

I  a  f  ^2t^i,x^in  ■*■  ^ltCi,x2m  +  ®3r^i,x^f  ■*■  I 
f  6^  •  dx 

""  36^(11,^  y  +  Ci,*  )Z|  +  S2rCi,yZi 


(62) 
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M  F 


12  2 

b-lr-li-l 


1‘  1 
dy  +  1 

1. 

1 

b  1 

dx  1 

.  J 

0 


(63) 


The  Ritz  polynomial  displacement  functions  do  not  satisfy  the 
boundary  conditions  of  a  plate,  singly  supported  on  four  sides. 

The  displacement  of  the  mid  surface  at  x=a  and  y=b  are  not  zero. 
The  rotations  and  at  x=0  and  y=0  are  equal  to  zero.  The  Ritz 

polynomial  represents  a  plate,  fixed  or  clarq>ed  on  two  sides  at 
x=0  and  y=0.  This  satisfies  the  conditions  that  ^=u=0  at  x=0  and 
^=v=0  at  y=0.  Also,  at  the  free  edge  x=a,  u  and  *P  do  not  equal 
zero  but  'P,x=0  and 'P,y=0  and  at  y=b,  ^,x=0  and^,y=0.  This  satisfies 
boundary  conditions  61  and  62.  Boundary  condition  63  is  satisfied 
by  the  condition  that  w=0  at  x=0  and  y=0.  Only  the  geometric 
boundary  conditions  are  satisfied  by  the  Ritz  Polynomial 
displacement  assun^tions. 

The  governing  equations  can  be  transformed  into  a  linear 
system  of  equations  as  follows: 


[»HBi  ^  -  (Oj  J  i  -  O..H^ 

[cj  Ip) 


(64) 


M  will  be  [2*n2+1]  x  [3*n2]  matrix  while  the  Ai,  Bi,  and  Ci  matrix 
will  be  3*n2  vector  and  the  right  hand  matrix  will  be  a  2*n2+i 
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vector.  The  x,  y  and  z  integrals  of  the  governing  equations  can 
be  evaluated  directly  leaving  a  system  of  equations  as  function  of 
the  Ai,  Bi  and  Ci  coefficients  only.  Guassian  elimination  can  be 
used  to  solve  for  the  coefficients  Ai,  Bi,  and  Ci.  The 
coefficients  are  then  used  to  generate  the  displacements,  strains 
and  stress  throughout  the  laminate. 

The  polynomial  order  N,  for  a  Ritz  polynomial  displacement, 
must  include  enough  terms  to  represent  the  complex  strain 
derivatives.  A  higher  order  of  N  is  required  for  a  Ritz 
polynomial  than  a  Sine/Cosine  function  to  properly  represent  the 
displacement.  An  N  of  six  creates  36th  order  polynomial  and  73  X 
108  [M]  matrix.  The  Ritz  polynomial  solution  does  not  converge 
monotonically  on  the  correct  solution  as  N  is  increased.  The 
solutions  generated  oscillate  between  very  good  and  very  poor. 


Simply  Supported  0/60  Composite  Plate  Example 

The  same  sample  plate  described  in  Table  2  and  3  is  used  for 
the  Ritz  polynomial  angle  ply  exan^le  accept  that  the  02  is  60 

degrees.  An  N  of  6  was  used  creating  36th  order  polynomials  in 
each  direction.  The  displacements,  strains  and  stresses  of  a 
simply  supported  con^site  plate  with  0/60  plies  are  presented  in 
figures  68-82.  The  x  and  y  dependence  has  been  factored  out  to 
give  the  through  the  thickness  distributions. 
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Normalized  z 

)/60  Laminated  Plate 


Normalized  z 


Figure  74.  Epsilon  (ey)  for  a  0/60  Laminated  Plate 


Normalized  z 


Figure  75.  Epsilon  (e,)  for  a  0/60  Laminated  Plate 
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Figure  76. 


Gamma  (Yj^)  for  a  0/60  Laminated  Plate 


^ - 1 - 1 - 1 - h 


0 

Normalized  z 


Figure  77.  Gamma  (Ys)  ^or  a  0/60  Laminated  Plate 


Normalized  z 


Figure  78.  Gamma  (Yyz)  for  a  0/60  Laminated  Plate 


Sigma 


Normalized  z 


Figure  79.  Sigma  (o^)  for  a  0/60  Laminated  Plate 


8] 


-0.5  0  0.5 

Normalized  z 

Figure  84.  Tau  (Xyjt)  for  a  0/60  Laminated  Plate 

Ritz  Polynaaial  0/60  Composite  Plate  Discussion 

As  can  be  seen  in  the  displacement  functions  in  figures  69 
and  70,  an  N  value  of  6  creating  a  36^  order  polynomial 
representation  of  the  displacement  function,  falls  short  of  the 
Sine/Cosine  function  5^  order  polynomial  modeling  of  the 
displacement.  The  outer  matrix  strata  exhibit  displacements  that 
are  not  consistent  with  the  loading  and  the  boundary  conditions. 
The  strain  functions  exhibit  somewhat  better  results  but  still  are 
not  consistent  with  classical  solutions.  An  additional  problem 
arises  with  the  neutral  axis  of  the  plate  in  the  Ritz  polynomial 
representation.  The  integral  of  the  stress  over  the  area  of  the 
plate  should  produce  the  stress  in  each  direction.  Figures  79  and 
80  show  the  neutral  axis  remaining  at  the  mid  plane.  The  integral 
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of  the  stress  generates  an  in  plane  force  which  should  not  be 
present.  The  Ritz  polynomial  also  contain  the  Scune  discontinuous 
transverse  shear  strain  and  illogical  transverse  stress  result 
found  in  the  Sine /Cosine  solution.  The  above  problems  in  the 
Ritz  polynomial  solution  make  it  a  guestionabe  method  for  solving 
the  Strata  model. 
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VII.  ConcluBiona  and  R0C<mmeadation8 


Conclusions 

The  purpose  of  the  strata  model  is  to  more  accurately  predict 
the  strain  and  stress  distributions  with  in  an  individual 
con^site  ply  and  throughout  the  thickness  of  a  composite 
laminate.  The  model  is  expected  to  be  most  useful  when  the  fiber 
and  matrix  materials  used  have  dissimilar  stiffness  properties  or 
high  temperature  brings  that  about.  For  the  simple  cases  run  in 
this  document,  the  strata  model  predicts  fiber  strata  normal 
stresses  5  times  as  high  as  those  experienced  by  the  matrix 
strata.  The  model  generates  transverse  normal  stresses  o,  on  the 

same  order  of  magnitude  as  the  in  plane  shear  stresses.  This 
suggests  that  transverse  normal  stresses  should  not  be  assumed 
negligible,  as  in  some  models.  The  strata  model  assumption  of  a 
zero  z  gradient  in  the  fiber  strata  places  the  fiber  strata  in  a 
state  of  plain  strain.  This  produces  discontinuous  and  Yyz 

shear  strains.  This  anomaly  of  the  strata  model  can  not  be 
removed  without  making  the  solution  procedure  vastly  more 
complicated.  From  the  plain  strain  condition  in  the  fiber  strata, 
the  model  predicts  significant  transverse  shear  stresses  in  the 
fiber  strata.  Most  models  deem  transverse  shear  as  negligible. 
Overall  the  model  predicts  a  significant  variation  of  the  stresses 
through  the  thickness  of  an  individual  ply.  This  is  a  sharp 
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constrast  to  the  constant  strain  and  stress  predicted  by  current 
methods . 

The  strata  model  is  based  on  an  assumed  displacement  of  a 
polynomial  in  z  and  solves  for  the  x  and  y  dependent 
displacements.  The  strata  model  can  be  solved  with  assumed  x  and 
y  displacement  functions  of  either,  a  Sine/Cosine  series,  or  a 
Ritz  polynomial  series. 

The  Sin/Cosine  series  solution  used  in  this  work  limits  the 
application  to  cos^sites  with  only  cross  plies,  in  order  to 
factor  out  the  x  and  y  dependence  instead  of  integrating  them  out 
of  the  differential  equations.  The  Sine/Cosine  series  solution 
requires  a  moderate  polynomial  order  to  model  the  x  and  y 
dependence,  but  requires  a  full  series  expansion  to  represent 
single  loading.  Factoring  out  the  x  and  y  dependence  allows  the 
system  of  equations  to  be  solved  without  losing  any  displacment 
mode  shsapes.  A  double  series  expansion  of  Sin/Cosine  functions 
could  be  used  to  model  composites  with  angle  plies  but  would 
require  a  large  number  of  Sine/Cosine  integrals  to  integrate  out 
the  X  and  y  dependence. 

A  Ritz  polynomial  series  requires  a  much  higher  polynomial 
order,  but  allows  the  analysis  of  con^sites  with  angle  plies. 

This  is  possible  because  the  x  and  y  dependence  can  be  integrated 
out  of  the  governing  equations  without  losing  any  displacement 
modes.  An  extremely  high  order  Ritz  polynomial  is  required  to 
model  a  sin?)le,  two  ply,  composite  plate  and  obtain  results 
consistent  with  the  loading  and  other  models.  The  Ritz  polynomial 
solution  does  not  conerge  monotonically  on  the  correct  solution  as 
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the  order  N  is  increased.  The  con^utational  time  becomes 
prohibitive  as  the  polynomial  order  grows.  Ritz  polynomial 
displacement  function  assuiq>tions  are  not  recommended  to  solve 
strata  model  problems. 

Areas  for  Further  Research 

This  thesis  only  investigated  singly  supported  plates  with 
distributed  pressures.  How  other  types  of  boundary  conditions  can 
be  modelled  with  the  strata  model  should  be  investigated  further. 

Sin/Cosine  terms  in  a  double  series  expansion  should  be 
investigated  as  a  displacement  assumptions  for  the  strata  model. 
This  type  of  displacement  assun^tion  has  the  potential  to  enable 
the  analysis  of  con^sites  with  cross  plies,  without  requiring  an 
inordinate  polynomial  order. 

The  effects  of  the  order  of  the  displacement  polynomial  were 
investigated  in  a  very  limited  sense  in  this  work.  More  complex 
loading  and  an  increased  number  of  plies  would  require  higher 
order  polynomials  to  represent  the  deformation  and  strain. 

Further  research  in  this  area  is  necesscuiy. 

One  of  the  reasons  the  strata  model  was  developed  is  to 
better  model  dynamic  loading.  Dynamic  terms  can  be  included  in 
the  derivation  of  the  governing  equations  and  a  larger  more 
complex  system  of  equations  esteiblished  and  solved. 

The  through  the  thickness  stress  distributions  generated  by 
the  strata  model  should  be  investigated  as  a  causes  of  known 
failure  mechanisms.  The  unique  stress  and  strain  distributions 
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generated  by  the  strata  model  could  be  used  to  explain  known 
failure  mechanisms  in  current  composites. 

One  of  the  primary  advantages  of  the  strata  model  is  it 
allows  matrix  dominated  areas  to  deform  with  some  independence 
from  the  stiff er  fiber  strata.  In  high  ten^rature  applications 
the  matrix  material  can  become  viscoelastic  and  its  properties 
degrade.  Further  investigation  is  required  on  the  effects  of 
degrading  the  matrix  properties. 
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Appendix  A:  Mathnatica  Routine  to  Generate  the  Strain  Energy 
Terms 


(♦This  Routine  generates  the  strain  Energy  matrices  *) 

Clear [all] 

(♦Enter  the  stiffness  Matrices  ♦) 

s={{  si,  s3,sl0,  s5,  0,  0}, 

{  s3,  s2,sll,  s6,  0,  0}, 

{sl0,sll,sl2,sl3,  0,  0}, 

{  s5,  s6,sl3,  s4,  0,  0}, 

{  0,  0,  0,  0,  s7,  s9}, 

{  0,  0,  0,  0,  s9,  s8}} 

k={{kl,k2,k2,  0,  0,  0}, 

{kl,kl,k2,  0,  0,  0}, 

{kl,k2,kl,  0,  0,  0}, 

{  0,  0,  0,k3,  0,  0}, 

{  0,  0,  0,  0,k3,  0}, 

{  0,  0,  0,  0,  0,k3}, 

(♦  Enter  the  strain  vectors  ♦) 

em  =  {{umx},  {vmy},  {wmz},  {(uny+vmx)},  { ( umz+vnnx ) } ,  {(vmz+wmy)}} 

ef  =  -{{ufx},  {vfy},  {wfz},  {(ufy+vfx)},  {(ufz+vnnx)},  {(vfz+wfy)}} 

emt  =Transpose[em] 
eft  =Transpose[ef ] 

(♦  Create  the  Strain  energy  Matricies  ♦) 
MatrixStrataStrainEnergy=( 1/2 ) ♦emt.k.em 
FiberStrataStrainEnergy  = ( 1 /2 ) ♦eft . s . ef 
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Appendix  B:  Mathmatica  Routine  to  Solve  a  Homogeneous  Plate 
Clear [all] 


Np=5 

j=Range[0,Np] 

nplies=2 
plythick= . 5 

LaminateHeight=plythick*nplies 
hf={. 3535, .3535, .3535} 

Do  [nrt:[ [f ] ]=plythick-hf [ [f ] ] , {f ,nplis+l}] 

Do[zf [ [f ] ]=( (nplies/2)-(f-l/2)*plythick, {f ,nplies+l} ] 

zm={.5, .42675, .07325, -.07325, -.42675, -.5} 
zntt={. 5,  .07325, -.4267} 
zinb={42675 ,  - .  07325 ,  - .  5} 

volf=0 

volnF=l-volf 

a=10 

b=10 

p=Pi/a 

q=Pi/b 

theta={0,0} 

Pl=l 


Ef=  {20000000,20000000} 
vf=  {.35, .35} 

Gf=  {Ef[[l]]/(2*(l+vf[[ll])),Ef[[231/{2*(l+vf[[2]]))} 

Ef=  {300000,300000} 
v£=  {.4, .4} 

Gf=  {Eni[[l]]/(2*(l+vm[[l]])),Em[[211/(2*(l+vm[[2]]))} 

zf 1 [ [ 1 ] ]=Laminateheight/2 

zfl[ [nplies+2] ]=-Laminateheight/2 

Be=plythick/ ( plythick-hf ) 

Do[al[m]=(  (zfl[  [m]  ]/ziitt[  [m]  ]  )-(zf  1[  [m+l ]  ]/zmb[  [m]  ] )  )* 

( (zmt[  [m]  ]*zinb[  [m]  ] ) /  (zinb[  [m]  ]-zmt[  [m]  ] ) ) ,  {m,nplies}  ] 

s=  Array[sf , {7,nplies}] 
k=  Array[km, {3,nplies+l} ] 


DO  [km[  1  ,in]=(Em[  [m]  ] / ( ( l-2*vm[  [m]  ] ) * ( l+vm[  [m]  ] ) ) * ( l-vm[  [m]  ] ) ) , 
{m,nplies+l}] 

DO  [lan[l,m]=(Em[[m]]/((l-2*vin[[m]])*(l+vm[[m]]))*(  vm[[m]])), 

{m,nplies+l}] 

DO  [kin[l,m]=(Em[  [m]]/(  (l-2*vm[  [m] ] )*(l+vm[ [m]] ) )*(l-2*vm[ [m]] )/2) , 
{m,nplies+l} ] 


91 


(♦Calculate  the  Fiber  strata  Engineering  Properties  *) 

DolEll[ [f ] ]=volf*Ef [ [f ] ]+volm*Em[ [f ] ] , {f ,nplies} ] 

Do[E33[  [f  ]  ]=volf*Ef  [  [f  ]  ]+volm*Ein[  [f  ]  ] ,  {f  ,nplies}  ] 

Do[E22[ [f ] ]=Em[ [f ] ]*Ef [ [f ] ]/{volf*Ef [ [f ] ]+volm*Em[ [f ] ] ) , 

{f ,nplies} ] 

Do[vl2[  [f  ]  ]=volf*vf  [  [f  ]  ]+volm*vin[  [f  ]  ] ,  {f  ,nplies}  ] 

Do[vl3[ [f ] l=volf*vf [ [f ] ]+volm*vm[ [f ] ] ,{f ,nplies} ] 
Do[v23[[f]]=vl2[[f]]*E22[[f]]/Ell[[f]],{f,nplies}] 

Do [ G1 3 [ [ f ] ] =vol f *Gf [ [ f ] ] +volm*Gm[ [ f ] ] , { f , npl ies } ] 

Do[G12[[f] ]=Em[[f ]]*Ef [[f ]]/(volf*Ef[(f ]]+volm*Em[ [f]]), 

{f ,nplies} ] 

Do[E23[[f]]=G12[[f ]],{f,nplies}] 

(♦Calculate  the  fiber  strata  stiffness  matrix  form  the  Eng  Con^) 


DO  [  V[[f]]=l- 

vl2[[f]]^(vl2[[f]]^E22[[f]]/Ell[[f]]+v23[[f]]^vl3[[f]])- 
vl3[ [f ] ]"2^1-v23[ [f ] ]^2^E33[[f ]]/E22[ [f ] ],{f,nplies}] 

DO  [sf[l,f]=(l-v23[[f]]-2^E33[[f]]/E22[[f]])^Ell[[f]]/V[[f]], 
{f ,nplies} ] 

DO  [sf[2,f]=(l-vl3[[f]]-2^E33[[f]]/Ell[[f]])^E22[[f]]/V[[fl], 
{f ,nplies} ] 

DO  [sf[3,f]=(vl2[[f]l+vl3[[f]]^E33[[fl]/E22[[f]])^ 
(E22[[f]]/V[[f]]),  {f,nplies}] 

DO  [sf[4,f]=(v23[If]]+vl2[[f]]^vl3[[f]]^E22[[f]]/Ell([f]])^ 
(E33[[f]]/V[[f]]),  {f,nplies}] 

Do  [sf[5,f]=G12[[f ]],  {f,nplies>] 

DO  [sf[6,f ]=G13[ [f ] ],  {f,nplies}] 

Do  [sf[7,f]=G12[[f ]],  {f,nplies>l 


(♦Calculate  the  rotated  fiber  strata  stiffness  matrix  (sr)  ♦) 


sr=Array [ sf r , { 1 3 , nplies } ] 


Do[sfr[l,f]  = 


Do[sfr[2,f]  = 


Do[sfr[3,f]  = 


Sin [ theta [ [f ] 
Cos [ theta [ [f ] 
Cos[theta[ [f ] 
Sin[theta[[f ] 
Cos[theta[ [f ] 
{f, nplies}] 
Sin[theta[ [f ] 
Cos[theta[ [f ] 
Cos[theta[ [f ] 
Sin [ theta [ [f ] 
Cos [theta] [f] 
{f, nplies} ] 
Sin] theta] [f] 
Cos ] theta] [f] 


]"2^(Sin]theta]]f ] ]"2^sf ]2,f]+ 
]"2^sf ]3,f ]+ 

]"2^(Cos]theta]]f  ]  ]''2^sf  ]l,f  ]+ 
]^2^sf ]3,f]+ 

]"2^Sin]theta]]f]]"2^sf]5,f], 

]''2^(Sin]theta]]f  ]  ]''2^sf]2,f  ]+ 
]''2^sf  [3,f]+ 

]"2^(Cos]theta]]f ]]"2^sf]l,f ]+ 
]"2^sf [3,f 1+ 

]"2^Sin]theta]]f]]"2^sf]5,f], 

]"2^(Sin]theta]]f ] ]"2^sf ]2,f ]+ 
]"2^sf ]3,f ]+ 
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Cos[theta[  [f  ]]'‘2*(Cos[thetaIIf  ]  ]"2*sf  [l,f  ]  + 
Sin[theta[[f ]]*2*sf[3,f]+ 

Cos [ theta[  [  f  ]  ] 2*Sin[ theta[  [  f  ]  ] '' 2*sf  [ 5 , f  ] , 

{f ,nplies} ] 

Dotsfr[4,f]  =  2*Cos[theta[[f]]]*Sin[theta[[f]]l* 

(2*Cos[theta[ [f ]] ]*Sin[theta[ [f ] ]]*sf[2,f]- 
Cos[theta[  [f  ]  ]  ]'‘2+Sin[theta[  [f  ]  ]  ]"2*sf  [5,f  ] , 

{f ,nplies}] 

Do[sfr[5,f]  =  2*Cos[theta[[f]]]*Sin[theta[[f]]]* 

(Sin[theta[ [f]]]"2*s£[2,f]+Cos[theta[ [f ] ] ]"2*sf [3,f ] )- 
2*Cos[theta[ [f ] ] ]*Sin[theta[ [£] ] ]* 
(Cos[theta[[f]]l"2*sftl,f]+Sin[theta[[f  ]]]''2*sf[3,f  ])  + 
Cos [ theta [ [ f ] ] ] *Sin [ theta [ [ f ] ] ] * 
(Cos[theta[[f]]]"2+Sin[theta[[f l]])*sf [3,f]], 

{f ,nplies}] 

Do[sfr[6,f]  =  -2*Cos[theta[[f]]]*Sin[theta[[f]]]* 

(Sin[theta[[f  ]]]^2*sf  [l,£]+Cos[theta[  [f  ]  ]  ]''2*sf  [3,f  ]  )- 
2*Cos[theta[ [f ] ] ]*Sin[theta[ [f ] ] ]* 
(Cos[theta[tf]]]"2*sf[2,f]+Sin[theta[[f ]]]^2*sf[3,f])+ 
Cos [ theta[ [ f ] ] ] *Sin[theta [ [ f ] ] ] * 
{Cos[theta[[f]]]"2+Sin[theta[[f ]l])*sf [5,f]], 

{f ,nplies}] 

Do[sfr[7,f ]  =Cos[theta[ [f ] ] ]"2*sf [6,f ]+Sin[theta[ [f ] ] ]"2*sf [7,f ] , 

{£/nplies} ] 

Dotsfr[8,f]  =Sin[theta[[f]]]"2*sf[6,f]+Costtheta[[f]]]"2*sf[7,f], 

{f ,nplies} ] 

Dolsfr[9,f]  =-(Cos[theta[[f]]]*Sin[theta[[f]]]*sf[6,f])+ 

Cos[theta[ [f ] ] ]*Sin[theta[ [f J ] ]*sf [7,f ], {f ,nplies} ] 

Do[sfrtlO,f]  =  Cos[theta[tf]]]^2*sf(3,f]+Sin[thetat[f]]]*sf[4,f] 

{f ,nplies>] 

Do[sfr[ll,f ]  =Sin[theta[ [f ]]]^2*sf[3,f]+Cos[theta[[f]]]"2*sf [4,f ], 

{£,nplies} ] 

Do[sfr[12,f]  =  sf [l,f],{f,nplies}] 

Do[sfr[13,f]  =  2*Cos[theta[[f]]l*Sin[theta[[f]l]*sf[3,f])+ 

2*Cos[theta[[f]]]*Sin[theta[[f]]]*sf[4,f], 

{f ,nplies}] 

(*  Put  the  stiffness  into  matrix  form  *) 

S=Table[SS[f ] , {f ,nplis} ] 

K=Table[KK[m] , {m,nplies+l} ] 

Do[SS[fl={{  sr[[l]][[fll,  sr[[311[[fl],sr[[10ll[[f]],sr[[5]],0,0}, 
{  sr[[3]][[f]],  sr[[2]][tf]],sr[[ll]][[f]],sr[[6]],0,0}, 
{sr[[10]][[f]],sr[[lll][[f]],sr[[12]][[f]],sr[[13]],0,0} 
{  sr[[5]]t[f]],  sr[[3]][[f]],sr[[10]][[f]],sr[[5]],0,0}, 
{  0,  0,0,0,  sr[[7]][[f]],sr[[9]][[f]]}, 

{  0,  0,0,0,  sr[[9]][[f]],sr[[8]][[f]]}} 

,{f,nplies}] 

Do[KK[m]={{  k[(l]][[m]],  k[[2]][Em]],k[[2]][[m]], 0,0,0}, 

{  k[[21l[[ml],  k[[l]][[mll,k[[2]l[[ml], 0,0,0}, 

{  k[[2]][[m]],  k[[2]][[m]],k[[l]][[m]], 0,0,0}, 

{  0,  0,  0,k[[3]][[m]],0,0}, 
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0,0,k[[3]][[m]],0}, 

0,0,0,k[[3]][[m]]}} 


{  0,  0, 

{  Of  Of 

,  {in,nplies+l}  ] 

(♦Enter  expressions  for  the  load  and  phi,  chi  and  w  functions*) 

P  =  Pl*Sin[p*xl*Sin[q*yl 

B1  =  Table[BB[i],{i,0,Np}] 

Cl  =  Table[CC[i],{i,0,Np}] 

phi  =  Bl*Cos[p*x]*Sin[q*y] 
chi  =  Cl*Sin[p*x]*Cos[q*y] 
w  =  A  *Sin[p*x]*Sin[q*y] 

(♦Generate  the  derivatives  of  phi,  chi  and  w  ♦) 

chix  =  D[chi,x] 
chiy  =  D[chi,y] 
chixx  =  D[chix,xl 
chiyy  =  D(chiy,y] 
chixy  =  D[chix,y] 

phix  =  D[phi,x] 
phiy  =  D[phi,y] 
phixx  =  D[phix,x] 
phiyy  =  D[phiy,y] 
phixy  =  D[phix,y] 

wx  =  D[w,x] 
wy  =  D[wy] 
wxx  =  D[wx,x] 
wyy  =  D[w,y] 
wxy  =  D[wx,y] 

(♦  Create  the  Differential  Equationr  for  a  Homogeneous  Plate*) 
edu  =  j  '  { 

Sum[ Integrate [  k[  [3,m]  ]*z^ ( j-1 ) , {z,znib[  [m]  ]  ,ziitt[  [m]  ]} ] , 

{m,nplies+l}]*  m/x  + 

Sum [ Integrate [sr[  [7,f  ]  ]*z''(  j-1) ,{z,zmt[  [f+l]  ]  ,zmb[  [f  ]  ]} ] , 

{f,nplies>]*  wx  + 

Sum[ Integrate [sr[ [9,f ] ]*z^( j-1) ,{z,zmt[ [f+1] ] ,zmb[ [f ] ]}] , 

{f,nplies}]*  wy) 

+  (Sum[-i*j* 

(Sum[ Integrate (  k[  [3,m]  ]*z''(  j+i-2),  {z,zmb[  [m]  ]  ,zmt[  [m]  ]>] , 

{m,nplies+l}]*  phi[[i+l]]+ 

Sum[lntegrate[sr[ I7,f ] ]*z^ ( j+i-2) , {z,zmt[ [f+1] ] ,zmb[ [f ] ]>] , 

{f,nplies}]*  phi[[i+l]]+ 

Sum[ Integrate [sr[ [9,f ] ]*z^( j+i-2) , {z,zmt[ [f+1] ] ,zmb[ [f ] ] } ] , 

{f,nplies}]*  chi[[i+l]]  )+ 

(Sum[ Integrate [sr[  [Iff]  ]*z'' ( j+i)  ,{z,zmt[  [f+1]  ] ,zmb[  [f  ]  ]} ] , 

{f ,nplies} ]+ 
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Suin[lntegrate[  k[  [l,m]  ]*z^  ( j+i) ,  {z,zinb[  [m]  ]  ,zint[  [m]  ] }  ] , 

{m,nplies+l} i*  phixx[ [i+1] ]+ 

2 *Sum[ Integrate [sr[  [5,£]  ]*z^  ( j+i) ,  {z,znrt: [  [f+1]  ]  ,znib[  [  f  ]  ] }  ] , 

{f,nplies}]*  phixy[ [i+1 ] ]+ 

(Sum[ Integrate [sr[  [4,f  ]  ]*z^(  j+i) ,  {z,zmt[  [f+1  ]  ]  ,zinb[  [ f  ]  ]}  ] , 

{f,nplies>]  + 

Sum[lntegrate[  k[  [3,m]  ]*z^ ( j+i) ,  {z,zinb[  [m]  ]  ,zmt[  [m]  ]} ] , 

{m,nplies+l} ] )*  phiyy[ [i+1] ]+ 

Suin[ Integrate [sr[  [5,f]  ]*z^(  j+i) ,  {z,zmt[  [f+1  ]  ]  ,zinb[  [£]]}], 

{f,nplies}]*  chixx[[i+l]]  + 

Sum[ Integrate [sr[  [6,f  ]  ]*z^  ( j+i) , {z,zmt[  [f+1  ]  ]  ,zinb[  [£]]}], 

{f,nplies}]*  chiyy[[i+l]]  + 

(Sum[lntegrate[sr[  [3,f]  ]*z''(  j+i)  ,{z,zmt[  [f+1]  ]  ,zmb[  [£]]}] , 

{f ,nplies} ]+ 

Sum[Integiate[sr[  [4,f  ]  ]*z''  ( j+i) , {z,zmt[  [f+1  ]  ]  ,zinb[  [f  ]  ]} ] , 

{f,nplies}]  + 

Suin[Integrate[  k[  [2,m]  ]*z'‘ ( j+i) ,  {z,zinb[  [m]  ]  ,znrt[  [m]  ]} ] , 

{m,nplies+l} ]+ 

Sum[ Integrate [  k[  [3,inj  ]*z^ ( j+i) ,  {z,zinb[  [m]  ]  ,zntt[  [m]  ]} ] , 

{m,nplies+l} ]  )*  chixy[[i+l]] 


,  {i,0,Np}]  ) 


edv  =  j*( 

Sum [ Integrate [k[ [3,m] ]*z^( j-1) , {z,zmb[ [m] ] ,zmt[ [m] ]}] , 

{m,nplies+l} ]*  wy  + 

Sum[ Integrate [sr[ [9,f ] ]*z^ ( j-1) , {z,zmt[ [f+1 ] ] ,zmb[ [f ] ]} ] , 

{f^nplies}]*  wx  + 

Sum[ Integrate [sr[ [8,f ] ]*z" ( j-1 ) , {z,zmt[ [f+1 ] ] ,zmb[ [f ] ]}] , 

{f,nplies)]*  wy) 

+  (Sum[-i*j* 

( Sum [ Integrate [  k[[3,m]  ]*z'' ( j+i-2) , {z,zmb[  [m]  ],zmt[  [m]]}], 

{m,nplies+l} ]*  chi[[i+l]]+ 

Sum[ lntegrate[sr[  [8,f  ]  ]*z''(  j+i-2) ,  {z^zmt[  [f+1  ]  ]  ,zmb[  [f  ]  ] }  ] , 

{f,nplies}]*  chi[[i+l]]+ 

Sum[lntegrate[sr[  [9,f  ]  ]*z''(  j+i-2) ,  {z,zmt[  [f+1]  ]  ,zmb[  [f  ]  ]}] , 

{f^nplies}]*  chi[[i+l]]  )+ 

(Sum[ Integrate [sr[ [2,f ] ]*z^ ( j+i) , {z,zmt[ [f+1] ] ,zmb[ [f ] ] } ] , 

{f ,nplies} ]+ 

Sum[lntegrate[  k[  [ l,m]  ]*z''  ( j+i) ,  {z,zmb[  [m]  ]  ,zmt[  [m]  ] } ] , 

{m,nplies+l> ]*  chiyy[ [i+1 ] ]+ 

2 *Sum[ Integrate [sr[ [6, f] ( j+i) , {z,zmt[ [f+1] ] ,zmb[ [ f ] ]} ] , 

{f,nplies}]*  chixy[ [i+1 ] ]+ 

(Sum[ Integrate [sr[  [4,f  ]  ]*z''(  j+i) ,  {z,zmt[  [f+1]  ]  ,zmb[  [f  ]  ]}  ] , 

{f,nplies}]  + 

Sum[Integrate[  k[ [3,m]  ]*z" ( j+i) , {z,zmb[ [m] ] ,zmt[ [m] ] } ] , 

{m,nplies+l}] )*  chixx[ [i+1] ]+ 

Sum [ Integrate [sr[ [6,f ] ]*z^ ( j+i) , {z,zmt[ [f+1 ] ] ,zmb[ [f ] ]} ] , 

{f,nplies}]*  phiyy[[i+l]]  + 

Sum[ Integrate [sr[  [5,f  ]  ]*z''  ( j+i) ,  {z,zmt[  [f+1  ]  ]  ,zmb[  [f  ]  ]} ] , 

{f,nplies}]*  phixx[[i+l]]  + 

(Sum[ Integrate [sr[ [3,f ] ]*z^( j+i) , {z,zmt[ [f+1] ],zmb[ [f ]]}] , 

{f ,nplies} ]+ 
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Sum [ Integrate [sr[  [4,f  ]  ]*z'' ( j+i) ,  {z,zmt[  [f+1  ]  ]  ,zinb[  [f  ]  ]}  ] , 

{f,nplies}]  + 

Sum [ Integrate [  k[ [2,m] ]*z^ ( j+i) , {z,zmb[ [m] ] ,zmt[ [m] ]} ] , 

{m,nplies+l} ]+ 

Sum[  Integrate [  k[  [3,m]  ]*z''  ( j+i) ,  {z,zmb[  [m]  ]  ,zmt[  [m]  ] }  ] , 

{m,nplies+l} ]  )*  phixyi[i+l]] 


,  {i,0,Np}]  ) 


edw  = 


-{ Sum [ Integrate [  k[ [ 3,m] ] , {z,zmb[ [m] ] ,zmt[ [m] ] } ] , {m,nplies+l} ]+ 
Sum[ Integrate [sr[ [7,f]] , {z,zmt[ [f+1] ] ,zmb[ [ f ] ] } ] , {f ,nplies} ] ) 

*  wxx  + 

2 *Sum[ Integrate [sr[ [ 9,f ] ] , {z,zmt[ [f+1 ] ] ,zmb[ [f ] ]} ] , {f ,nplies} ] 

*  WTT/ 

( Sum [ Integrate [  k[ [ 3,m] ] , {z,zmb[ [m] ] ,zmt[ [m] ] } ] , {m,nplies+l} ]+ 
Sum[ Int‘^grate[sr[  [8,f  ]  ] ,  {z,zmt[  [f+1  ]  ]  ,zmb[  [f  ]  ] }  ] ,  {f  ,nplies} ] ) 

wyy  + 


Sum[-i( 

Sum[ Integrate I  k[ [3,m]  ]*z^(i-l),{z,zmb[ [m] ] ,zmt[  [m] ]}] , 

{m,nplies+l}]*  wy  + 

Sum[ Integrate [sr[ [9,f ] ]*z^ ( j-1) , {z,zmt[ [f+1 ] ] ,zmb[ [f ] ]} ] , 

{f,nplies}]  *  (phix[ [i+1 1 ]+chiy[ [i+1 ] 

Sum[ Integrate [sr[ [7,f ] l*z^{ j-1 ) , {z,zmt[ [f+1 ] ] ,zmb[ [f ’ ] } ] , 

{f,nplies}]  *  phix[[i+l]]  + 

Sum[Integrate[sr[  [8,f  ]  ]*z''  ( j-1 ) ,  {z,zmt[  [f+1  ]  ]  ,zmb[  [f  ]  ] }  ] , 

{f/nplies}]  *  chiy[[i+l]] 


,{i,0,Np}] 


])  + 


(*  Construct  a  system  of  equations  for  the  Diff  Equations  *) 

m  =  Array[mm,{2*(Np+l)+l,2*(Np+l)+l>] 

f  edu=edu  /Cos  [  p*  x  ]  *  S  in  [  q*y  ] 

Do[M[[j+l,i+l]]  = 

Coefficient [fedu[ [ j+1 ] ] ,BB[i] ] , { jrOfNp} , {i,0,Np} ] 

Do[M[ [ j+l,i+l+Np+l] ]  = 

Coefficient[fedu[ [ j+1] j,CC[i]],{j,0,Np},{i,0,Np}] 

Do[M[ [ j+l,2*(Np+l)+l] ]  = 

Coefficient [fedu[ [ j+1] ] /A] j,0,Np}, {i,0,Np} ] 

f edv=edu/Cos [ q*x ] *Sin [ p*y ] 

Do[M[[ j+Np+2,i+l]]  = 

Coefficient [fedu[ [j+1 ] ] ,BB[i] ] , { j,0,Np}, {i,0,Np} ] 

Do[M[ [ j+Np+2,i+l+Np+l] ]  = 

Coef ficient[fedu[ [ j+1 ] ] ,CC[i] ] / { j,0,Np}, {i,0,Np} ] 

DO[M[[ j+Np+2,2*(Np+l)+l]]  = 

Coef f icient [ f edu [ [ j+1 ] ] , A] , { j , 0 , Np} , { i , 0 , Np} ] 

f edw=edu/Sin [ p*x ] *Sin [ q*y ] 

Do[M[ [2*(Np+l)+l,i+l] ]  = 

Coef ficient [ fedu[ [ j+1 ]] ,BB[ i] ],{ j , 0,Np} , {i, 0 ,Np} ] 

Do[M[ [2*(Np+l)+l,i+l+Np+l] ]  = 

Coefficient [fedu[ [ j+1] ] ,CC[i] ] , { j,0,Np}, {i,0,Np} ] 
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DO[M[ [2*(Np+l)+l,2*(Np+l)+l] ]  = 

Cc5efficient[fedu[  [  j+1]  ]  ,A] ,{  j,0,Np>,  {i,0,Np}  ] 

(*  Create  the  answer  vector  (bm)  for  the  ra*sol=bm  system  *) 

bm  =  Array  [bbm, {2*(Np+l),l}] 

DO  [bm[[i]]=0,  {i,2*(Np+l)}] 

Do  [bm[ t2*(Np+l)+l] ]=P1,  {i,2*(Np+l)}] 

(*  Solve  the  system  of  Equations  for  the  sol  vector  *) 

sol=LinearSolve [ m, bm ] 

(*  Determine  the  Bl,  Cl,  A  coefficients  *) 

Do[BB[i]=sol[ [i+1] , {i,0,Np} ] 

Do[CC[i]=sol[[i+Np+2},{i,0,Np}] 

A=sol[ [2*(Np+l)+l] ] 

(*  Creat  the  Displacement  functions  for  the  fiber  strata  *) 

Uf [z_]=Table[uff [z] [f ] ,{f ,nplies}] 

Vf  [z_i=Ted)le[vff  [z]  [f  ]  ,{f  ,nplies}] 

Wf [z_]=Table[wff [z] [f ] , {f ,nplies}] 

Do[uff[z_][f]  =  Sum[zf[[f]]"i*BB[i],{i,0,Np},{f,nplies>] 
Do[vff[z_] [f ]  =  Sum[zf[[f]]^i*CC(i],{i,0,Np},{f,nplies}] 
Do[wff[z_][f]  =  A,{f,nplies}] 

(*  Create  the  Displacement  Functions  for  the  matrix  strata  *) 

Um[z_]=Table[umm[z] [m] , {m,nplies+l> ] 
vm[z_i=Table[vmm[zi [m] , {m,nplies+l}] 

Wm[  z_]  =Table  [  wmm[  z  ]  [  m] ,  {m,  nplies+1  >  ] 

Do[umm[z_][f]  =  Sum[(al[[m]]+Bem[[m]]*z)"i*BB[i], 

{i,0,Np>, {m, nplies+1 > ] 

Do[vmm[z_] [f ]  =  Sum[ (al[ [m] ]+Bem[ [m]  ]*z)^i*BB[i] , 

{i,0,Np},{m, nplies+1 } ] 

Do[wmm[z_] [f ]  =  A, {m, nplies+1 } ] 

(*  Create  the  Strain  Functions  for  the  matrix  strata  *) 

efx[z_]=Table[effx[z] [f ] , {f ,nplies>] 
efy[z_]=Table[effy[z] [f ] , {f ,nplies>] 
efz[z_]=Table[effz[z] [f ] , {f ,nplies}] 
efxy[z_]=Table[effxy[z] [f ] , {f ,nplies} ] 
efxz[z_]=Table[effxz[z] [f ] , {f ,nplies}] 

ef yz[z_] “Table [ eff yz[z] [£] r {f ^nplies}] 

Do  [effx[z_][f l=p*Uf[zl[[f ]l,{f ,nplies}] 

DO  [effy[z_] [f l=p*Vf [z] [ [ f ] ] , {f ,nplies> ] 

Do  [effz[z_] [f ]=0, {f ,nplies} ] 

DO  [effx[z_][f]=q*Uf[z][[f]]+p*Vf[z3[[f]],{f,nplies>] 

DO  [effy[z_][f]=D[Uf[z][[f]],z]+p*Wf[z][[f]],{f,nplies}] 

DO  [effy[z_][f]=D[Vf[z]I[f]],z]+q*Wf[z][[f]],{f,nplies>] 
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(*  Create  the  Strain  Functions  for  the  matrix  strata  *) 


emx[z_l=Tcd3le[emmx[2] [m] , {m,nplies+l} ] 
eiiiy[z_]=Table[eimny[z]  [m] ,  {m,nplies+l>  ] 
emzf z_]=Table[emmz[zj [m] , {m,nplies+l} ] 
emxy  [  z_]  =Tal>le  [  emmxy  [  z  ]  [m] ,  {m,  nplies+1 }  ] 
emxz [ z_] “Table [emmxz [ z ] [m] , {m, nplies+1 } ] 
eayz [ z_] “Table [ emmyz [ z ] [ “ ] # {®/ nplies+1 } ] 

Do  [emffix[z_] [f ]“-p*Umt2] [ [f ] ] ,{f ,nplies}] 

Do  [ emmy [ z_] [ f ]“-q*Vm[ z ] [ [ f ] ] / { f , nplies} ] 

Do  [emmz[z_][f ]“0,{f,nplies}] 

Do  [eiranx[z_][f l“q*Um[zl[ [m]]+p*Vf tz][ [m] ],{m,nplies+l}] 

Do  [emny  [ z_]  [ f ]“D tUm[ z ] [ [m] ] , z ]+p*Wf [ z ] [ [m] ] , {m, nplies+1 > ] 

Do  [eimry[z_]  [f  l=D[Vm[z]  [  [m]  ]  ,z]+q*wf  [z]  [  [m]  ] , {m,nplies+l} ] 

(*  Create  the  Strain  vector  for  the  fiber  strata  *) 

Fstrain  [  z_]  =TeU3le  [  f  strain  [z][f],{f,  nplies  }  ] 

Do[fstrain[z_] [f ]“{efx[z] ,efy[2] [f ] ,efz[z] [f ] »efxy[z] [f ] # 

efxz[zl[f],efyz[z]if]},{f, nplies}] 

(*  Create  the  Strain  vector  for  the  matrix  strata  *) 

Mstrain [ z_] =Table [ mstrain [ z ] [ m] , {m, nplies+1 } ] 

Do[mstrain[z_]  [f  ]={emx[z]  [m]  ,emy[z]  [m]  ,emz[z]  [m]  ,emxy[2]  [m] , 

emxz  [  z  ]  [m] ,  enyz  [  z  ]  [m] } ,  {m,  nplies+1 }  ] 

(*  Calculate  the  Stress  in  each  stratum  by  Stiff ness*strain  *) 

Fstress [ z_] “  Tcible [ f stress [z][f],{f, nplies } ] 

Do[ f stress [ z_] [f]“  S[ [f ] ] .Fstrain[z] [ [f ] ] , {f, nplies} ] 

Mstress[z_]“  Table [mstress[z] [m],{m, nplies+1}] 

Do[mstress[z_] [m]=  S[ [m] ] .Fstrain[z] [ [m] ], {m, nplies+1 } ] 


(*  Place  the  stress  con^x^nentes  into  plotable  form  *) 

smx[2_]“Table[smxxtz]  [m] ,  {m,  nplies+1}] 
smy  [  z_]“TeQ)le[  smyy  [  z  ]  [m] ,  {m,  nplies+1}  ] 
smz[z_]“Table[smzz[z] [m] , {m,nplies+l}] 
smxy [ z_] “Table [ smxyy[ z ] [m] , {m, nplies+1 } ] 
smxz [ z_ ] “Table [ smxzz [ z ] [m] , {m, nplies+1 } ] 
smyz [ z_] “Table [ smyzz [ z ] [ m ] , {m, nplies+1 } ] 

sfx[z_]“Ted3le[sfxx[z]  [f  ] ,  {f, nplies}] 
sfy[z_]“TeU&le[sf3^[z]  [f  ] ,  {f, nplies}] 
sfz[z_]“Table[sfzz[z] [f], {f, nplies}] 
sf xy [ 2_]“Table [ sfxyy [ z ] I f ],{ f , nplies } ] 
sfxz[z_]=Table[sfxzz[z] [f ] , {f, nplies} ] 
sf yz ( z_] “Table [ sf yzz [ z ] [ f ] , { f , nplies} ] 
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Dotsinxx[z_]  [m]  =Mstress[z]  [  [m]  ][[!)],  {m,nplies+l>] 
Do[sntyy[z_i  [m]  =Mstress[z]  [  [m]  ]  [  [2]  ] ,  {m,nplies+l}] 
Do  [  smzz  [  z_]  [  m]  =Mstress  [z][[in]][[3]]/{m,  nplies+1 }  ] 
Dotsinxyy[z_]  [m]=Mstress[z]  [  [m]  ]  [  [4]  ] ,  {m,nplies+l}  ] 
Do[  smxzz  [  z_]  [m]=Mstress  [z][[®]l[l5]],  {m,  nplies+1  >  ] 
Do[snyzz[z_]  [ini=Mstress[z]  [  [m]  ]  [  [6]  ],  {m, nplies+1}] 

Do[sfxx[z_] [f ]  =Fstress[z] [ [f ] ] [ [1] ], {f,nplies}] 
Do[sfyy[z_] [f ]  =Fstress[z] [ [f ] ] [ [2] ], {f ,nplies}] 
Do[sfzz[z_] [f ]  =Fstress[z] [[f]l[[3l] , {f ,nplies>] 
Do[sfxyy[z_] [f ]=Fstress[z] [ [f ] ] [ [4] ], {f ,nplies}] 

Do[sfxzz[z_] [f )=Fstress[z] [ [f ] ] [ [5] ] , {f ,nplies>] 

Do[sfyzz[z_] [f ]=Fstress[z] [ [f ] ] [ [6] ] , {f ,nplies}] 
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Appendix  C:  Mathmatica  Routine  for  Strata  Model  Problems 
Clear [all] 

Np=5 

j=Range[0,Np] 

nplies=2 
plythick= . 5 

LaininateHeight=plythick*nplies 
hf={. 3535,. 3535} 

Do  [ntt[  [f  ]  ]=plythick-hf  [  [f  ]  ] ,  {f  ,nplis+l} ] 

Do[zf [ [f ] ]=( (nplies/2)-(f-l/2)*plythick, {f ,nplies+3 } ] 

ZIB={.5,  .42675,. 07325, -.07325, -.42675, -.5} 
zmt={. 5, .07325, -.4267} 
zmb={42675,-.07325,-.5} 

volf=l. 4142/2 

voln^l-volf 

a=10 

b=10 

p=Pi/a 

q=Pi/b 

theta={0,0} 

Pl=l 

Ef=  {20000000,20000000} 
vf=  {.35,. 35} 

Gf=  {Ef[[l]]/(2*(l+vf[[l]])),Ef[[2]]/(2*(l+vf[[2]]))} 

Ef=  {300000,300000} 
vf=  {.4, .4} 

Gf=  {Em[[lll/(2*(l+vm[[l]l)),EmC[2]]/(2*(l+vm[[2ll))} 

zf  1  [  [  1  ]  ]=LeUiiinateheight/2 

zf 1 [ [nplies+2 ] ]=-Lamlnateheight/2 

Be=plythick/ ( plythick-hf ) 

Do[al[m]=( (zfl[ [m] ]/zmt[ [m] ] )-(zfl[ [m+l] ]/zmb[ [m] ] ) )* 

( (zmt[ [m] ]*zmb[ [m] ] ) / (zmb[ [m] ]-zmt[ [m] ] ) ) , {m,nplies} ] 

s=  Array[sf , {7,nplies}] 
k=  Array [ km, { 3 , nplies+1 } ] 

DO  [km[l,m]=(Em[[m]]/((l-2*vm([m]])*(l+vm[[m]]))*(l-vm[[m]])), 

{m, nplies+1}] 

DO  [km[l,m]=(Em[[m]]/((l-2*vm[[m]])*(l+vm[[m]]))*(  vm[[m]])), 

{m, nplies+1 }  ] 

DO  [km[l,m]=(Em[[m]]/((l-2*vm[[m]])*(l+vm[[m]]))*(l-2*vm[[m]])/2), 
{m,  nplies+1 } ] 
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(♦Calculate  the  Fiber  strata  Engineering  Properties  *) 


Do[El 1 [ [ f ] ]=volf *Ef [ [ f ] ] +volm*Em[  [ f ] ] , {f ,nplies> ] 

Do[E33[  [f  ]  ]=volf*Ef  [  [f  ]  ]+volin*Em[  (f  ]  ]  ,{f  ,nplies}  ] 
Do[E22[[f]]=Em[[f ]]*Ef[[f]]/(volf*Ef[[f]]+volm*Em[[f]]), 

{f ,nplies}] 

Do[vl2[ [f ] ]=volf*vf [ [f ] ]+volm*vm[ [f 1 ] , {f ,nplies}] 

Do[vl3[  [f  ]  ]=volf*vf  [  [f  ]  ]+volni*vm[  [f  ]  ] ,  {f  ,nplies}] 
Do[v23[[f]]=vl2[[f]]*E22[[f]]/Ell([f]],{f,nplies>] 

Do[G13 [ [ f ] ]=volf *Gf [ [ f ] ] +volm*Gm[ [ f ] ] , {f , nplies} ] 

Do[G12[[f ]]=Em[[f ]]*Ef[[f]]/(vol£*Ef( [f ] ]+volm*Em[ [f ] ] ) , 

{f, nplies}] 

Do[E23[[f]]=G12[ [£]],{£, nplies}] 

(♦Calculate  the  fiber  strata  stiffness  matrix  form  the  Eng  Con^) 
Do  [  V[[f]]=l- 

vl2[[f]]^(vl2[[fl]^E22[[f]]/Ell[[fll+v23[[f]]^vl3[[f]l)- 
vl3[[f]]*2^1-v23[[f]]''2^E33[[f]]/E22[[f]],{f, nplies}] 

DO  [sf[l,f]=(l-v23[[f]]"2^E33[[f]]/E22[[f]])^Ell[[f]]/V[[f]], 

{f, nplies} ] 

Do  [sf[2,f]=(l-vl3([f]]"2^E33[[f]]/Ell[[f]])^E22[[f]]/V[[f]], 

{f, nplies}] 

DO  [sf [3, f]=(vl2[[f]]+vl3[[f]]^E33[[f]]/E22[ [£]])♦ 
(E22[[f]]/V[[f]]),  {f, nplies}] 

Do  [sf [4, f]=(v23[[f]]+vl2[[f]]^vl3([f]]^E22[(f]]/Ell[ [£]])♦ 
(E33[[f]]/V[[f]]),  {f, nplies}] 

Do  [sf[5,f]=G12[[f]],  {f, nplies}] 

Do  [sf[6,f]=G13[[f]],  {f, nplies}] 

Do  [sf[7,f]=G12[[f ]],  {f, nplies}] 

(♦Calculate  the  rotated  fiber  strata  stiffness  matrix  (sr)  ♦) 
sr=Array [ sf r , {13 , nplies} ] 

Do[sfr[l,f]  =  Sin[theta[[f  ]]''2^(Sin[theta[[fl]''2^sf[2,f]+ 
Cos(theta[[f ]]"2^sf [3,f ]+ 

Cos [ theta[ [ f ] ] ^2^ (Cos [ theta[ [ f ] ] "2^sf [ 1 , f ]+ 
Sin[theta[[f  ]l''2*sf  [3,f  ]+ 

Cos ( theta[  [  f  ]  ] ''2*Sin[ theta[  [ f  ]  ] '‘2*sf  [ 5 ,  f  ] , 

{f, nplies}] 

Do[sfr[2,f]  =  Sin[theta((f]]"2^(Sin(theta[(f]]''2*sf[2,f]+ 

Cos[theta[[f ]]"2^sf [3,f ]+ 

Cos [ theta[  [  f  ]  ] ''2* (Cosi theta[  ( f  ]  ] ''2^sf  [  1 , f  ]+ 
Sin[theta[ [f ] ]"2*sf [3,f ]+ 

Cos  [  theta  [  [  f  ]  ]  ^  2  *Sin  [  theta  [  [  f  ]  ] '' 2  *s  f  [  5 ,  f  ] , 

{f, nplies}] 

Do[sfr[3,f]  *  Sin[theta([f]]''2^(Sin[theta[[f]]^2^sf[2,f]+ 

Cos[theta[[f ]]^2*sf [3,f ]+ 


Cos[theta[[f]]''2*(Cos[theta[[f  ]  ]''2*sf[l,f  ]+ 
Sin[theta[[f ]]^2*sf[3,f]+ 

Cos[theta[[f ]]]"2*Sin[theta[[f ] ]"2*sf [5,f ], 

{f ,nplies} ] 

Do[sfr[4,f]  =  2*Cos[theta[[£]]]*Sin[theta[[f]]]* 

(2*Cos[theta[[£]]]*Sin[theta[[f]]]*sf[2,f]- 
Cos [ theta[ [fill *2+Sin[ theta[ [ f 1 1 r2*sf [ 5 , f 1 , 
{£,nplies} 1 

Do[sfr[5,fl  =  2*Cos[theta[[flll*Sin[theta[[flll* 

(Sin[theta[[f  lir2*sf  [2,f  l+Cos[theta[[f  lir2*sf[3,fl)- 
2*Cos [theta [[fill *Sin[ theta[ [fill* 

(Cos[theta[  [f  1 1  r2*sf  [l,f  1+Sin[theta[  [f  1 1  r2*sf  [3,f  1 )+ 
Cos[theta[ [f 1 1 l*Sin[theta[ [f 1 1 1* 
(Cos[theta[[fllK2+Sin[theta[[f  ll])*sf[3,fll, 

{f ,nplies}l 

Do[sfr[6,fl  =  -2*Cos[theta[[fl]l*Sin[theta[[fll]* 

(Sin[theta[[flir2*sf[l,fl+Cos[theta[  [f  1 1  K2*sf  [3,f  1  )- 
2*Cos [ theta [[fill *Sin[ theta[ [fill* 

(Cos[theta[[f  lir2*sf  [2,f  1+Sin[theta[[f  lir2*sf[3,f  1)  + 
Cos [ theta [ [ f l 1 1 *Sin [ theta [[fill* 

(Cos [ theta[  [ f  1 1 1 ''2+Sin[ theta [  [ f  1 1 1 ) *sf  [5, f  1 1 , 

{f ,nplles} 1 

Do[sfr[7,f  1  =Cos[theta[  [flir2*sf[6,fl+Sin[theta[  [f  lir2*sf[7,f  1, 

{f ,nplies} 1 

Do[sfr[8,f 1  =Sin[theta[ [flll"2*sf[6,f l+Cos[theta[ [f 1 1 r2*sf [7,f 1 , 

{f ,nplies}l 

Do[sfr[9,fl  =-(Cos[theta[[f]ll*Sin[theta[[f lll*sf[6,f])+ 

Cos[theta[[fl]l*Sin[theta[[f]ll*sf[7,fl,{f,nplies}l 

Do[sfr[10,fl  *  Cos[theta[[flir2*sf[3,f ]+Sin[theta[[flll*sf[4,f] 

{f ,nplies}l 

Do[sfr[ll,fl  =Sin[theta[[flir2*sf(3,fl+Cos[theta[[flir2*sf[4,fl, 

{f ,nplies}l 

Do[sfr[12,fl  =  sf [l,f l,{f ,nplies}l 

Do[sfr[13,fl  =  2*Cos[theta[[flll*Sin[theta[[flll*sf[3,fl)+ 

2*Cos[theta[ [f 1 1 l*Sin[theta[ [f 1 1 l*sf [4,f 1 , 

{f ,nplies} 1 

(*  Put  the  stiffness  into  matrix  form  *) 

S=Table[SS[f l,{f,nplis}l 

K=Table[KK[m] , {m,nplies+l} 1 

Do[SS[fl={{  sr[[lll[[fll,  sr[[3ll[[fll,sr[[10]l[[fll,sr[[511,0,0>, 
{  sr[[3ll[[fll,  sr[[2ll[[fll,sr[[llll[[fll,sr[[6ll,0,0}, 
{sr[ [101 1 [ [f 1 1, sr[ [111 l[[f]l,sr[ [121 l[[fll,sr[ [1311,0,0} 
{  sr[[5ll[[fll,  sr[[3ll[[fll,sr[[10ll[[fll,sr[[5ll,0,0}, 
{  0,  0,0,0,  sr[[711[[fll,sr[[9ll[[fll}, 

{  0,  0,0,0,  sr[[911[[fll,sr[[8ll[[fll>} 

,{f,nplies}l 

Do[KK[ml={{  k[[lll[[mll,  k[[21l[[ml],k[[2ll[[mll,0,0,0>, 

{  k[[2ll[[mll,  k[[lll[[mll,k[[2ll[[mll, 0,0,0}, 

{  k[[2]l[[mll,  k[[2ll[[m]l,k[[lll[[mll, 0,0,0}, 

{  0,  0,  0,k[[3ll[[mll,0,0}, 
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0,0,k[[3]][[m]],0}, 

0,0,0,k[[3]][[ml]}} 


{  0,  0, 

{  0,  0, 

,{m,nplies+l>] 

(♦Enter  expressions  for  the  load  and  phi,  chi  and  w  functions*) 

P  =  Pl*Sin[p*xl*Sin[q*yl 

B1  =  Table[BB[i],{i,0,Np}] 

Cl  =  Table[CC[i],{i,0,Np}] 

phi  =  Bl*Cos[p*x]*Sin[q*y] 
chi  =  Cl*Sin[p*x]*Cos[q*y] 
w  =  A  *Sin[p*x]*Sin[q*y] 

(♦Generate  the  derivatives  of  phi,  chi  and  w  ♦) 

chix  =  D[chi,x] 
chiy  =  D[chi,y] 
chixx  =  D[chix,x] 
chiyy  =  D[chiy,y] 
chixy  =  D[chix,yj 

phix  =  D[phi,x] 
phiy  =  D(phi,y] 
phixx  =  D[phix,x] 
phiyy  =  D[phiy,y] 
phixy  =  D[phix,y] 

wx  =  D[w,x] 
wy  =  D[wy] 
wxx  =  D[wx,x] 
wyy  =  D[w,y] 
wxy  =  D[wx,y] 

(♦  Create  the  Differential  Equations  for  a  Coiqx>site  Plate  ♦) 
edu  =  j*Be( 

Sui&[ Integrate [  k[  (3,m]  ]*(al[  [m]  ]+bem[  [m]  )*z)''  ( j-1) , 

{z,zinb[  [m]  ]  ,zint[  [m] ]}] , {m,nplies+l> ]  )*wx  + 

+  (Sum[-i*j* 

(Sum[ Integrate (  k[ (3,m] ]*(al( [m]]+bem[ [m] ]*z)'‘ ( j+i-2), 

{z,zinb[  [m]  ]  ,zntt[  [m]  ]}] ,  {m,nplies+l}]*  phi[  [i+1]  ]+ 

(Sum[sr[[l,f]]*hf[[f]]*zf[[fJ]-(j+i),{f,nplies>]+ 

Sum[  Integrate [  k[  [  1  ,m]  ]  ♦  ( al [  [m]  ]+bem[  [m]  ]  ♦z ) ''  ( j+i) , 

{z,ziiib[  [m]  ],zmt[  [m]  J}],{ m,nplies+l}]) ♦phixx [  [i+1]  ]+ 
2*Sum[sr[[5,f]]*hf[[fl]*zf[[f]]^(j+i),{f,nplies}]*phixy[[i+l]]+ 
(Suin[sr[l4,f]]*hf[[f]]*zf[[fJ]''(j+i),{f,nplies}]  + 

Sum [ Integrate [  k[  [3,m]  ]*(al[[m]  ]+bem[  [m]  ]*z)''(  j+i) , 

{z,zmb[ [m] ],zmt[ [m] ]>], {m,nplies+l}] ) ♦phiyy [ [i+1] ]+ 
Sum[sr[[5,f]]*hf[  [f]]*zf[[f]]''(  j+i),{f,nplies>]*chixx[[i+l]]+ 
Sum[sr[ [6,f ] ]*hf [ [f ] ]*zf [ [f ] ]" ( j+i) , {f ,nplies}] ♦chiyy [ [i+1] ]+ 
(Sum[sr[[3,f]]*hf[[f]]*zf[[f]]Mj+i)r{frnplies}]+ 
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Sum[sr[[4,f]]*hf[[f]]*zf[[fl]"( j+i),{f,nplies}]  + 

Sum[ Integrate [  k[  [2,ni]  ]*(al[  [m]  ]+bem[  [m]  ]*2)"  (:+i) , 
{z,zmb[[m]],zmt[[in]l>l,{m,nplies+l}]+ 

Suin[  Integrate [  k[  [  3  ,m]  ]  *  ( al [  [m]  ]+bem[  [m]  ]  *z  r  (  d+i ) ,  , , .  ^ , 

{ z ,  zinb[  [ m]  ] ,  zmt  [  [  m]  ] }  ] ,  {m,  nplies+1 }  ]  )  *chxxy  [  [  i+l  ]  ] 

,  {i,0,Np}l  ) 


©d.v  "  3  ^  ( 

Sumf lntegrate[k[ [3,in] ]*(al[ [m] ]+bem[ [m] ]*z)  ( j-1) / 

{z,znib[[m]],ziitt[[m]]>],{m,  nplies+1}]*  wy  + 

+  (Suni[-i*j*  .  .  _ 

(Sum[ Integrate [  k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*zr  (3+1-2) , 

{z,znib[  [m]  ]  ,zittt[  [m]  ]}]  ,{m,nplies+l}]*chi[  [i+l]  ]  + 

Sum[sr[[8,fll*h£[[f]l*zf[[fir(j+i-2),{f,nplies>]*chi[[i+l]]+ 

Sum[sr[[9,f]]*hf[[f]]*zf[[£]]"{j+i-2),{f,nplies>] 

*cni[ [1+1 ] ] )+ 

(Suni[sr[[2,£]l*h£[[£ll*z£[[£l]"(j+i),{f»nplies}]+ 

Sum[ Integrate [  k[ [l^m] ]*(al[ [m] ]+bem[ [m] ]*z)  (3+i). 

{z,znib[  [m]  ]  ,zmt[  [m]  ] } ] ,  {m,nplies+l}  ]*  chiyy[  [i+l]  ]  + 

2*Sum[sr[[6,£]]*h£[[£]]*z£[[£]]"(j+i),{f/nplies}l*chixy[[i+l]]+ 
(Suin[sr[[4,£]l*h£[[£]]*z£[[£]]"(j+i),{f/nplies}]  + 

Suin[lntegrate[  k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*zr  (l+D  r  ^ 

{z,znib[  [m]  ] , zmt [  [m]  ] >  ] ,  {m,nplies+l}  ] ) *chixxt  [i+l ]  ]  + 
Suin[  sr  [  [  6 ,  £ ]  ]  *h£ [  [  £ ]  ]  *z£ [  [  £ ]  ]  M  j+i )  /  { f  / nplies }  ]  *Phiyy [  I i+J  1 1  + 

Sum[sr[[5,£]l*h£[(£]]*z£[[£]]"(3+i),{ffnplies>]*phixx[[i+l]]+ 
(Sum[sr[ [3, £]]*h£[[£]]*z£[t£]]"(j+i),{f, nplies} ]+ 

Suin[ sr[  [4, £]  ]*h£[  (£]]*z£[  [£]  ]"(3+i)  nplies}]  + 

Suin[ Integrate [  k[  [2,m]  ]*(al[  [m]  ]+bemt  [m]  ]*z)  (3+i)r 
{z,zinb[  [m]  1  ,zmt[  [m]  ] }  ] ,  {m,  nplies+1 }  ]+ 

Suin[ integrate [  k[  (3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*zr  (D+i)  r  r  r  , 
{z,zinb[[m]],znit[[m]]}],{m,nplies+l}]  )*phixy[  [i+l]  ] 

,  {i,0,Np}]  ) 


eclw  — 

-( Sum [ Integrate [k[ [3/m] ] , {z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ]+ 
Sum[sr[I7,£]]*h£[[£]],{£, nplies}])  *  wxx  + 

2*Sum[sr[[9,£]]*h£[[£]],{£, nplies}]  .  „  ^  t 

( Sum [ Integrate [k[ [3,m] ] , {z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l}]+ 
Sum[sr[[8,£]]*h£[[£]],{£, nplies}])  *  wyy  + 

Sum[-i( 

Sum[  lntegrate[ k[  [  3 ,m]  ] *  ( al [  [m]  j  +bem[  [ m]  ] *z ) ''  ( i-l ) , 

{z,zmb[lm]],zmt[[ml]}],{m,nplies+l}]*  wy  + 

r {if  Of Np} ] 

(*  Construct  a  system  o£  equations  £or  the  Di££  Equations  *) 

m  =  Array[mm, {2*(Np+l )+lf 2*(Np+l )+l} ] 


£edu=edu/Cos [ p*x ] *Sin [ q*y ] 

CTO££icient[£edu[ [ j+1] ] f BB[i] ] f { jf 0,Np} , {if 0,Np} ] 

Do[M[ [i+lf i+Np+2] ]  =  .  „  „  , 

Coe££icient[£edu[f[  j+1  ]]fCC[i]Jf  {  j,0,Np} ,  {if  0,Np}  ] 
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Do[M[[j+l,2*(Np+l)+l]]  = 

Coefficient [fedu[ [ j+1] ] ,A] , { j, 0,Np} ,{i,0,Np} ] 

f edv=edu/Cos [ q*x ] *Sin [ p*y ] 

Do[M[[ j+Np+2,i+l]]  = 

Coefficient [fedu[ [ j+1] ] ,BB[i] ] , { j,0,Np}, {i,0,Np} ] 
Do[M[[ j+Np+2,i+Np+2]]  = 

Coefficient [ f edu [ [ j+1 ]]#CC[i]],{jrO, Np} ,{1,0, Np} ] 
Do[M[[ j+Np+2,2*(Np+l)+l]]  = 

Coefficient [ f edu [ [ j+1 ])/A],{j, 0 ,Np> , {i , 0 ,Np} ] 

f edw=edu/Sin [ p*x ] *Sin [ q*y ] 

Do[M[[2*(Np+l)+l,i+l]]  = 

Coefficient [f edu [ [ j+1] ] ,BB[i] ] , { j/0,Np>, {i,0,Np>] 
Do[M[[2*(Np+l)+l,i+Np+2]]  = 

Coefficient [ fedu [ [ j+1 ] ] , CC[ i] ] , { j , 0,Np} , {i, 0 ,Np} ] 
Do[M[[2*(Np+l)+l,2*(Np+l)+ll]  = 

Coefficient [ fedu [ [ j+1] ] /A], { j,0,Np>, {i,0,Np} ] 

(*  Create  the  answer  vector  (bm)  for  the  m*sol=bm  system  *) 

bm  =  Array  [bbm, {2*(Np+l) , 1} ] 

Do  [bm[[i]]=0,  {i,2*(Np+l)}] 

Do  [bm[[2*(Np+l)+l]]=Pl,  {i,2*(Np+l)}] 

(*  Solve  the  system  of  Equations  for  the  sol  vector  *) 

sol=LinearSol ve [ m, bm ] 

(*  Determine  the  Bl,  Cl,  A  coefficients  *) 

Do[BB[i]=sol[[i+l],{i,0,Np>] 

Do[CC[i]=sol[ [i+Np+2},{i,0,Np}] 

A=sol[[2*(Np+-l)+l]] 


(*  Creat  the  Displacement  functions  for  the  fiber  strata  *) 

Uf [z_]=Table[uff [z] [f ] , {f ,nplies}] 

Vf  [2_]=Teible[vff  [z]  [f  ]  ,{f  ,nplies}] 

Wf  [z_]=Tcd)le[wff  [z]  [f  ]  ,{f  ,nplies}] 

Do[uff[z_][f]  =  Sum[zf[[f]]"i*BB[i],{i,0,Np},{f,nplies>] 
Do[vff[z_] [f ]  =  Sum[zf [ [f] ]^i*CC[i],{i,0,Np},{f,nplies>] 
Do[wff [z_] [f ]  =  A, {f ,nplies} ] 

(*  Create  the  Displacement  Functions  for  the  matrix  strata  *) 

Umtz_]=Ted)le[umm[z]  [m]  ,{m,nplies+l}] 

Vm[  z_]  =Table  [  vmm[  z  ]  [  m] ,  {m,  nplies+1 }  ] 

Wm[  z_]  “Table  [  wiiim[  z  ]  [  m] ,  {m,  nplies+1 )  J 
Do[umm[z_][f]  =  Sum[(al[  [m]]+Bem[lm]]*z)"i*BB[i], 

{ i , 0 , Np} , {m, nplies+1 } ] 

Do[vmm[z_][f  ]  =  Sum[(al[[m]]+Bem[lm]]*z)'‘i*BB[i], 

{ i , 0 , Np} , {m, nplies+1 } ] 
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Do[wnm[2_]  [f  ]  =  A,  {in,nplies+l}] 

(*  Create  the  Strain  Functions  for  the  matrix  strata  *) 

efx[z_]=Table[effx[z] [f ] , {£,nplies}] 
efy[z_]=Table[effy[zi [f ] , {£,nplies}] 
e£z[z_]=Table[e££z[z] [£] , {£,nplies}] 
e£xy[z_]=Table[e££xy[z] [£] ,{£,nplies} ] 
e£xz[z_]=Table[ef£xz[z] [f ] r {£/nplies}] 
e£yz[z_i=Table[e££yz[z] [£] , {£,nplies> ] 

DO  [e££x[z_] [£]=-p*U£[z] [ [£] ] , {£,nplies}] 

DO  [e££y[z_][£]=-q*V£[z][[£]],{£,nplies}] 

Do  [e££ztz_] [f ]=0# /Hplies}] 

Do  [e££x[z_][£]=q*U£[z][[£]]+p*V£[z][[£]],{£,nplies}] 

Do  Ie£fy[z_][£]=D[Uf[z]Il£]],z]+p*W£[z][[£]],{f,nplies>] 

Do  [e££y[z_] [£]=D[V£[z] [ [£] ] ,z]+q*W£[z] [ [£] ] , {£,nplies} ] 

(*  Create  the  Strain  Functions  for  the  matrix  strata  *) 

emx [ z_] =Table [ emmx [ z ] [ m] , {m, nplies+1 } 3 
emy[z_]=Table[emmy[z] [m] , {m,nplies+l}] 
emz[z_]=Teible[emmz[z]  [m] ,  {m,nplies+l}] 
emxy [ z_] =Table [ emmxy [ z 3 [ m 1 , {m, nplies+1 } ] 
emxz[z_3=Table[emmxz[z3 [mj , {m,nplies+l} 3 
emyz [ 2_3  =Table [ emmyz [ z  3 [mj , {m, nplies+1 } ] 

Do  [emmx[z_3 [£3=-p*Um[z3 [[£33 , {£,nplies}3 
Do  [emmy[z_3[£3=-q*Vm[z3[[£33,{£,npliesH 

Do  [emmz[z_3 [f 3=0, {£ ,nplies}3 

Do  [ emmx [ z_3 [ £  j  =q*Um[ z  3 [ [ m3  3  +p*V£ [ z  3 [ [m3  3 , {m, nplies+1 } 3 
Do  [emray[z_3[£3=D[Um[z3[[m33,z3+p*W£[z3[[m33,{m, nplies+1} 3 
Do  [eminy[z_3[£3=D[Vta[z3[[m33,z3+q*W£[z3[[m33,{m,nplies+l}3 

(*  Create  the  Strain  vector  for  the  fiber  strata  *) 

Fstrain[z_3=Table[fstrain[z3 [£3 ,{£ ,nplies>3 
Do[£strain[z_3 [£3={e£x[z3 ,e£y(z3 [£3 ,e£z[z3 [£3 ,e£xy[z3 [£3 , 

e£xz[z3[£3,e£yz[z3[f3>,{£,nplies}3 

(*  Create  the  Strain  vector  for  the  matrix  strata  *) 

Mstrain[z_3=Table[mstrain[z3 [m3 , {m, nplies+1} 3 
Do[mstrain[z_3  [£3={emx[z3  [mj  ,eniy[z}  [m3,emz[z3  [m3  ,emxy[z3  [m3 , 

emxz [ z 3 [m3 , emyz [ z 3 [ m3 } , {m, nplies+1 } 3 

(*  Calculate  the  Stress  in  each  stratum  by  Sti££ness*strain  *) 

Fstress[z_3=  Table[£stress[z3 [£3,{£,nplies}3 
Do[£stress[z_3 [£3=  S[ [£3  3 .Fstrain[z3 [ [£3  3 , {£,nplies}3 

Mstress[z_3=  Table [ms tress [z} [ra},{m, nplies+1} 3 
Do[mstress[z_3 [m3=  S[ [m3  3 .Fstrain[z3 [ [m3  3 , {m,nplies+l} 3 
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(*  Place  the  stress  con^nentes  into  plotable  form  *) 

smx[2_]=Teible[sinxx[z]  [m] ,  {m,nplies+l}] 
smy[z_]=Table[snyy[z] [m] , {m,nplies+l}] 
smz [ z_]=Table [ smzz [ z ] [m] , {m, nplies+1 } ] 
smxy [ z_] =Table [ smxyy [ z ] [ m] , {m, nplies+1 } ] 
smxz[z_]=Table[smxzz[z] [m] , {m, nplies+1 } ] 
smyz[z_]=Table[siiiyzz[z]  [m] ,  {m,  nplies+1 }  ] 

sf X [ z_] =Table [ s  f xx [ z ] [ f ] , { f , nplies } ] 
sfy[z_]=Table[sfyy[z] [f ] , {f ,nplies}] 
sfz[z_]=Table[sfzz[z][f],{f, nplies}] 
s  f xy [ z_ ] =Table [ s  f xyy [ z ] [ f ] , { f , nplies } ] 
sfxz[z_]=Table[sfxzz[z][£],{f, nplies}] 
sfyz [ z_]=TableI sfyzz [ z ] [ f ] , {£, nplies} ] 

Do[smxx[z_] [m]  =Mstress[z] [ [m] ] [ [1] ] , {m, nplies+1}] 

Do[ smyy [ z_] [m]  =Mstress [ z ] [ [m] ] [ [2 ] ] , {m, nplies+1} ] 

Do[ smzz [ z_] [m]  =Mstress [ z ] [ [m] ] [ [ 3 ] ] , {m, nplies+1} ] 

Do[ smxyy [ z_] [m]=Mstress [ z ] [ [m] ] [ [ 4 ] ] , {m, nplies+1} ] 

Do t  smxzz [ z_] [m]=Mstress [ z ] [ [m] ] [ [ 5 ] ] , {m, nplies+1 } ] 

DO [ srayzz [ z_] [ m] =Mstress [ z ][[ m] ][[ 6 ]], {m, nplies+1 } ] 

Do[s£xx[z_][f]  =Fstress[z] [[£]][[!]],{£, nplies}] 
Do[sfyy[z_][£]  =Fstress[z] [[£]][ [2] ], {f, nplies} ] 

Do[ sf zz [ z_] [ £ ]  =Fstress [z][[f]][[3]],{f, nplies} ] 
Do[sfxyy[z_] [f]aFstress[z] [[£]][ [4] ],{£, nplies}] 
Dots£xzz[z_][£]=Fstress[z] [[£]][ [5] ],{f, nplies}] 
Do[sfyzztz_][f]=Fstress[z][[f]][t6]],{f, nplies}] 
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Appendix  D:  Mathmatica  Routine  with  Ritz  Polynomial 
Clear [all] 


Np=5 

j=Range[0,Np''2  ] 

nplies=2 
plythick= . 5 

LaminateHeight=ply thick*  npl ies 
hf={. 3535, .3535} 

Do  [mt[[f]]=plythick-hf[ [fl],{f,nplis+l}] 

Do[zf [ [f ] ]=( (nplies/2)-(f-l/2)*plythick,{f ,nplies+l}] 

zmF{.5, .42675, .07325, -.07325, -.42675, -.5} 
zmt={. 5, .07325, -.4267} 
zmb={42675,-.07325,-.5} 

volf=l. 4142/2 

volm=l-volf 

a=10 

b=10 

p=Pi/a 

q=Pi/b 

theta={0,0} 

Pl=l 


zf  1  [  [  1  ]  ]=Lcuninateheight/2 

zfl[  [nplies+2]  ]=-Leuiiinateheight/2 


Be=plythick/ ( plythick-hf ) 

Do[al[m]=((zfl[[m]]/zmt[ [m] ] )-(zfl[ [m+1] ]/zmb[ [m] ] ) )* 

( (zmt[  [m]  ]*zinb[  [m]  ] ) / (zmbt  [m]  ]-zntt [  [m]  ] ) ) ,  {m,nplies}  ] 

Ef=  {20000000,20000000} 
vf=  {.35,. 35} 

Gf=  {Ef[[l]]/(2*(l+vf[(l]])),Ef[[2]]/(2*(l+vf[[2]]))} 

Ef=  {300000,300000} 
vf=  {.4,. 4} 

Gf=  {Em[[l]]/(2*(l+vm[[ll])),Em[{2]3/(2*(l+vm[[2]]))} 


s=  Array[sf , {7,nplies} ] 
k=  Array [ km, { 3 , nplies+1 } ] 


Do 

Do 

Do 


[km[ l,m]=(Em[  [m]  ]/( ( l-2*vm[  [m]  ]  )*(l+vm[  [m]  ] )  )*(l-vm[  [m]  ] ) ) , 

{m, nplies+1}] 

[km[l,m]=(Em[  [m] ]/( (l-2*vm[  [mll)*(l+vm[  [m]  ] )  )*(  vm[  [m]  ] ) ), 

{m, nplies+1}] 

[km[l,m]=(Em[[m]]/((l-2*vm[  [m]  ]  )*(l+vm[  [m]  ] )  )*(l-2*vm[  [m]  ]  )/2) , 
{m, nplies+1}] 
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(♦Calculate  the  Fiber  strata  Engineering  Properties  *) 


Do[Ell[ [f ] ]=volf*Ef [ [f ] ]+volm*Em[ (f ] ] , {f ,nplies} ] 

Do[E33[  [f  ]  ]=volf*Ef  [  [f  ]  ]+volm*Ein[  (f  ]  ] ,  {f  ,nplies}  ] 

Do[E22[ [f ] ]=Em[ [f ] ]*Ef [ [f ] ]/(volf*Ef [ [f ] ]+volm*Em[ [f ] ] ) , 
{f,nplies}] 

Do[vl2[ [f ] ]=volf*vf [ [f ] ]+volm*vm[ [f ] ] , {f ,nplies} ] 

Do[vl3[  [f  ]  ]=volf*vf  [  [f  ]  ]+volm*vin[  (f  ]  ] ,  {f  ,nplies}  ] 
Do[v23[[fl]=vl2([f]]*E22[(f]]/Ell[[f]l,{f,nplies>] 

Do[G13( [ f ] ]=volf*Gf [ [ f ] l+volm*Gm[ [f ] ] , {f ,nplies} ] 
Do[G12[[f]]=Ein[[f  ll*Ef[[fl]/(volf*Ef([f]]+volin*Em[[f  ]]), 

{f ,nplies} ] 

Do(E23[[f ]]=G12[[f l],{f,nplies}] 

(♦Calculate  the  fiber  strata  stiffness  matrix  form  the  Eng  Con*) 


Do  [  V[[f]l=l- 

vl2[(f]]^(vl2[[f]]^E22[(f]]/Ell[[f]]+v23[[f]]^vl3[[f]])- 
vl3[[f l]*2^1-v23[[f ]]"2^E33[[f ]]/E22[[f ]],{f ,nplies}] 


Do  [sf[l,f]=(l-v23[[f]]"2^E33[[f]]/E22I[f]])^Ell[[f]]/V[[f]], 
{f ,nplies} ] 

DO  [sf[2,f]=(l-vl3[(f]]"2^E33[[fl]/EllI[f]])^E22[[f]]/V[[fl], 
{f ,nplies}] 

DO  [sf[3,fl=(vl2[[f]l+vl3([f]]^E33[[f]]/E22[[f]])^ 
(E22[[f]]/V[(f]]),  {f,nplies}] 

DO  [sf (4, f]=(v23[[f]]+vl2[[f]]^vl3[[f]]^E22[[f]]/Ell[ [£]])♦ 
(E33((£]]/V[[f]]),  {f,nplies>] 

Do  [s£[5,f]=G12([f ]],  {f,nplies>] 

Do  [sf(6,fl=G13[[f 11,  {f,nplies}l 
Do  [sf[7,fl=G12[[£ll,  {f,nplies}l 


(♦Calculate  the  rotated  fiber  strata  stiffness  matrix  (sr)  ♦) 


sr=Array [ sf r , { 1 3 , nplies  > 1 


Do[sfr[l,fl  = 


Do(sfr[2,fl  = 


Do[sfr[3,fl  * 


Sin[theta[ [f 1 
Cos[theta[(f 1 
Cos[theta[ [f 1 
Sin[theta[ [f 1 
Cos[theta[[f 1 
{f, nplies} 1 
Sin [ theta [ [f 1 
Cos[theta[ [f 1 
Cos[theta[[fl 
Sin[theta[ [f 1 
Cos[theta[ [f 1 
{f, nplies}] 
Sin[theta[ [f 1 
Cos[theta[[f] 


r2*(Sin[theta[[f  ir2^sf  [2,f  1  + 

r2^sf  t3,fl  + 

r2^(Cos[theta[[f ll*2^sf [l,f 1+ 

r2*sfl3,fl+ 

r2*Sin[theta[(fir2^sf[5,f], 

r2*(Sinltheta[  [f  1  r2^sf  [2,f  1+ 
r2*sf  (3,f  1  + 

K2*(CosIthetal[f]r2^sf[l,f  1+ 
r2*sf[3,fl+ 

r2^Sin(theta([f  ir2*sf[5,f  1, 

r2*(Sin(theta[[f  1  r2*sf(2,f  1+ 
r2*sf  [3,f  1+ 
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Cos[theta[[f]]^2*(Cos[theta[[f]]*2*sf[l,f]+ 
Sin[theta[[f ] ]*2*sf [3,f]+ 

Cos[theta[ [f ] ] ]"2*Sin[theta[ [f ] ]^2*sf [5,f ] , 

{f ,nplies} ] 

Do[sfr[4,f]  =  2*Cos[theta[[f]]]*Sin[theta[[f]]]* 

(2*Cos[theta[ [f ] ] ]*Sin[theta[ [f ] ] ]*sf [2,f ]- 
Cos[theta[  [f  ]  ]  ]''2+Sinttheta[  [f  ]  ]  ]^2*sf  [5,f  ] , 

{f ,nplies} ] 

Do[sfr[5,f]  =  2*Cos[theta[[f]]]*Sin[theta[[f]]]* 

(Sin[theta[[f ] ]]"2*sf[2,£]+Cos[theta[ [f ] ] ]"2*sf[3,f ] )- 
2  *Cos [ theta [ [ f ] ] ] *Sin [ theta [ [ f ] ] ] * 

(Cos[theta[[f ]]]"2*sf[l,£]+Sin[theta[ [f ] ] ]"2*sf [3,f ] )+ 
Cos [ theta[ [ f ] ] ] *Sin[ theta [ [ f ] ] ] * 
(Cos[theta[[f]]]"2+Sin[theta[[f]]])*sf[3,f]], 

{f ,nplies} ] 

Do[sfr[6,f]  =  -2*Cos[theta[[f]]l*Sin[theta[[f] ]]* 

(Sin[theta[[f ]] ]"2*sf [l,f ]+Cos[thata[ [f ]]]"2*sf [3,f ] )- 
2*Cos[theta[ [f ] ] ]*Sin[theta[ [f ] ] ]* 

(Cos[theta[  [f  ]  ]  ]"2*sf  [2,f  ]+Sin[theta[  [f  ]  ]  ] ''2*sf  [3,  f  ] )  + 
Cos [ theta [ [ f ] ] ] *Sin ( theta [ [ f ] ] ] * 

(Cos[theta[ [f ] ] ] ^2+Sin[theta[ [f ] ] ] )*sf [5,f ] ] , 

{f ,nplies} ] 

Do[sfr[7,£]  =Cos[theta[ [f ]]]"2*sf[6,£]+Sin[theta[ [f ]] ]"2*sf[7,f ], 

{f ,nplies} ] 

Do[sfr[8,f ]  =Sin[theta[ [f]]]^2*sf[6,f]+Cos[theta[ [f ]]]"2*sf[7,f ], 

{f ,nplies}] 

Do[sfri9,f]  =-(Cos[theta[[£]]]*Sin(theta[[f]] J*sf[6,f])+ 

Costtheta[ [f ] ] ]*Sin[theta[ [f ] ] ]*sf [7,f ] , {f ,nplies} ] 

Do[s£r[10,f]  =  Cos[theta[[f]]]"2*sff3,f]+Sinttheta[[f]]]*sf [4,f] 

{f ,nplies) ] 

Do[sfr[ll,f  1  =Sin[theta[  [f  ]  ]  ]^2*sf  [3,f  ]+Cos[theta[  [f  ]  ]  ]'‘2*sf  [4,f  ] , 

{f ,nplies} ] 

Do[sfr[12,f]  =  sf [l,f],{f,nplies}] 

Do[sfr[13,f]  =  2*Cos[theta[[f]]]*Sin[theta[[f]]]*sf[3,f])+ 
2*Cos[theta[[f]]]*Sin[theta[[f]]]*sf[4,f], 

{f ,nplies} ] 

(*  Put  the  stiffness  into  matrix  form  *) 

S=Table [ SS [ f ] , { f , nplis } ] 

K=Table[KK[m] , {m,nplies+l} ] 

Do[SS[f]={{  sr[[l]][[f]],  sr[[3]l[[f]],sr[[10]][[f]],sr[[5]],0,0}, 
{  sr[[3]][[f]],  sr[[2]][[f]],sr[[ll]][[f]],sr[[6]],0,0}, 
{sr[[10]][[f]],sr[[ll]][[fl],sr[[12]][[f]],sr[[13]],0,0} 
{  sr[[5]][[f]],  sr[[3]][[f]3,sr[I10]][[f]],srt[5]],0,0}, 
{  0,  0,0,0,  sr[[7]][[fl],sr[[9]]t[f]]}, 

{  0,  0,0,0,  sr[[9]][[f]],sr[[8]][[f]]}} 

,{f,nplies}] 

Do[KK[m]={<  k[[l]][[m]],  k[[2]][[m]],k[[2]][[m]], 0,0,0}, 

{  k[[2]3[[m]],  k[[l]][[m]],k[I2]][[m]], 0,0,0}, 

{  k[[2]][[m]],  k[[2ll[[m]],k[[l]][[m]], 0,0,0}, 

{  0,  0,  0,k[[3]][tm]],0,0}, 
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{  0,  0,  0,0,k[l3]][[m]],0}, 

{  0,  0,  0,0,0,k[[3]][[m]]}} 

f {m,nplies+l} ] 

(♦Enter  expressions  for  the  load  and  phi,  chi  and  w  functions*) 

P  =  PI 

A1  =  Array[AAl,{Np,Np}] 

B1  =  Array[BBl,{Np,Np}] 

Cl  =  Array[CCl,{Np,Np}] 

A2  =  Table[AA[i],{i,Np"2}] 

B2  =  Table[BB[i],{i,Np''2}] 

C2  =  Table[CC[i],{i,Np"2}] 

DO  [Al[[n,l]]=A2[[l+(n-l)*Np]],{n,Np},{l,Np}] 

Do  [Bll[n,l]]=B2[[l+(n-l)*Np]],{n,Np},{l,Np}] 

Do  [Cl[[n,l]]=C2[[l+(n-l)*Np]],{n,Np},{l,Np}] 

phi  =  Sum  [Sum[Bl[[l,n]]*x''n*y"l,{n,NP}],{l,Np}] 
chi  =  Sum  [Sum[Cl[[l,n]]*x"n*y"l,{n,NP}],{l,Np}] 
w  =  Sum  [Sum[Al[[l,n]]*x"n*y"l,{n,NP}],{l,Np}] 

(♦Generate  the  derivatives  of  phi,  chi  and  w  ♦) 

Do  [chixl[i]  =  D[chi[[i]],x],  {i,Np^2}] 

DO  [chiylli]  =  D[chi[[i]],y],  {i,Np"2}] 

Do  [chixxl[i]  =  D[D[chi[[i]],x],x],  {i,Np"2}] 

DO  [chiyyl[i]  =  D[D[chi[ [i] ] ,y] ,yl ,  {i,Np"2>] 

Do  [chixyl[i]  =  D[D[chi[ [i] ] ,x] ,y] ,  {i,Np"2}] 

Do  [phixl[i]  =  D(phi[[i]],x],  {i,Np^2}] 

DO  [phiyl[i]  =  D[phi[[l]],y],  {i,Np-2}] 

DO  [phixxl[i]  =  D[D[phi[[i]l,x],x],  {i,Np-2}] 

Do  [phiyyl[il  =  D[D[phi[[i]],y],y],  {i,Np''2}] 

Do  [phixyl[i]  =  D[D[phi[[i]],x],yl,  {i,Np-2}l 

DO  [wxl[i]  =  D[w[[i]],x],  {i,Np-2>] 

DO  [wyl[i]  =  D[w[Ii]),y],  {i,Np-2}) 

DO  [wxxl[i]  =  D[D[w[[i]],x],x],  {i,Np"2}] 

DO  [wyyl[i]  =  D[D[w[[i]],y],y],  {i,Np"2}] 

Do  [wxyl[i]  =  D[D[w[[i]],x],y],  {i,Np^2}] 

(♦  Create  the  Differential  Equations  for  a  Con^site  Plate  ♦) 

edu  =  Integrate [ 

(Sum[ 

i*Be( 

Sum[ Integrate [  k[  [3,m]  ]*(al(  [m]  ]+bem[  [m]  ]*2) ''  (i-1 ) , 

{z,zmb[[m]],zmt[[m]]}],{m,nplies+l).]  )*wx[[i]]  + 

-i*  j* 

(Sum[ Integrate [  k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)''  ( j+i-2) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} 1 , {m,nplies+l} ]♦  phi[ [i] ]+ 
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(Sum[sr[[l,f]]*hf[[f]]*zf[[fl]Mj+i)r{f/nplies>]+ 

Sum [ Integrate [  k[ [1/m] ]*(al[ [m] ]+bem[ [m] ]*z) ^ ( j+i) , 

{z,zinb[  [m]  ]  ,zmt[  [m]  ]}  ] ,  {m,nplies+l} ]  )*phixx[  [i]  ]+ 
2*Sum[sr[[5,f]]*hf[[f]]*zf[[fll"( j+i),{f,nplies}]*phixyt[ill+ 
(Sum[sr[[4,f]]*hf[[f]]*zf[[f]]"(j+i),{f,nplies}]  + 

Sum [ Integrate [  k[ [3,m] ]*(al[[m] l+bem[ [m]  ]*z)^( j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ] ) *phiyy[ [i] ]+ 
Sum[sr[ [5,f ] ]*hf [ [f ] ]*zf [ [f ] ]"( j+i) ,{f ,nplies}]*chixx[ [i] ]+ 
Sum[sr[[6,f]]*hf[[f]]*zf[[f ]]"( j+i),{f,nplies}]*chiyy[[i]]+ 
(Sum[sr[[3,f]]*hf[[f]]*zf[[f]]"(j+i),{f,nplies}]+ 
SumLsr[[4,f]]*hfttf]]*zf[[f]]"(j+i),{f,nplies}]  + 

Sum[ Integrate [  k[ [2,m] ]*(al[ [m] ]+bem[ [m] ]*z)^ ( j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ]}] , {m,nplies+l} ]+ 

Sum[ Integrate [  k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)  '^  ( j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ]  )*chixy[ [i] ] 
,  {i,Np-2}l  ),  {x,0,a},{y,0,b}] 


edv  =  Integrate [ 

(Sum[ 

i*( 

Sum[ Integrate [ k [ [ 3 , m] ] * ( al [ [ m] ] +bem[ [ m] ] *z ) " ( i-1 ) , 

{z,zmb[ [m] ] ,zmt[ [m] 1}] ,{m,nplies+l}]*  wy[[i]]  + 

-i*  j* 

( Sum [ Integrate [  k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)''  ( j+i-2) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ]*chi[ [i+1 ] ]+ 
Sum[sr[[8,f]]*hf[[f]]*zf[[f]]"(j+i-2),{f,nplies}]*chi[[i]]+ 
Sumtsrt[9,f]]*hf[[f]]*zf[[f]]"(j+i-2),{f,nplies}] 

*chi[li]])+ 

(Sum[sr[[2,f]]*hf[[f]]*zf[[f]]-(j+i),{f,nplies}]+ 

Sum[ Integrate [  k[  [l,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)'‘(  j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ]*  chiyy[ [i] ]+ 
2*Sum[sr[[6,f]]*hf[[f]]*zf[[fJ]"(j+i),{f,nplies}]*chixy[[i]]+ 
(Sum[sr[[4,fl]*hf[[fl]*zf[[f3]"(j+i),{f,nplies}]  + 

Sum [ Integrate [  k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)''(  j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ] )*chixx[ [i] ]+ 
Sum[sr[[6,f]]*h£[:f]]*zf[[f]]"( j+i),{f,nplies}]*phiyy[[i]]+ 
Sum[sr[[5,f]]*hf[[f]]*zf[[f]]-( j+i),{f,nplies}]*phixx[[i]]+ 
(Sum[sr[[3,fl]*hft[f]]*zft[fll"(j+i),{f,nplies}l+ 

Sum[sr[[4,f  ]]*hf[  [f]]*zf  [[£]]''(  j+i),{f,nplies}]  + 

Sum[ Integrate [  k[  [2,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)''(  j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ]} ] , {m,nplies+l} ]+ 

Sum[  Integrate [  k[  [  3,m]  ]  *  (al [  [m]  ]+bem[  [m]  ]  *z )  '^  ( j+i) , 

{z,zmb[ [m] ] ,zmt[ [m] ] } 3 , {m,nplies+l} 3  )*phixy[ [±3  3 
,  {i^Np  2}3  )f  {X/ 0,a} , {y^ 0,b} ] 

edw  =  -Integrate [ 

( Sum [ Integrate [k[ [3,m3  3 , {z,zmb[ (m3  3 ,zmt[ [m3  3} 3 , {m,nplies+l} 3+ 
Sum[sr[[7,f33*hf[[f33,{f,nplies}3)  *  wxx[[i33  + 

2*Sum[sr[[9,f33*hf[[f33/{frnplies}3  *  wxy[[i33 

( Sum [ Integrate [k[ [3, m3  3 , {z^zmb[ [m3  3 ,zmt[ [m3  3 } 3 , {m,nplies+l} 3  + 
Sum[8r[[8,f33*hf[[f33,{f,nplies}3)  *  >.»yy[[i]3  + 

Sum[i( 
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Suin[ Integrate [ k[  [3,m]  ]*(al[  [m]  ]+bem[  [m]  ]*z)*  (i-1), 

{z,zn45[  [m]  ]  ,zmt[[m]  ]}]  ,{m,nplies+l}]*  wy  + 
,{i,Np-2}],  {x,0,a},{y,0,b}] 

(*  Construct  a  system  of  equations  for  the  Diff  Equations  *) 
m  =  Array  [  mm,  {2*Np ''2+1, 3*Np^2}  ] 
f edu=Expand [ edu ] 

Do[M[[j,i]]  =  CoefficientIfedu[[ j]],BB[i]],{j,Np-2},{i,Np"2}] 
Do[M[[j,i+Np-2]]  = 

Coefficient[fedu[ [ j]],CC[i]],{ j,Np"2},{i,Np"2}] 
Do[M[[ j,i+2*Np"2]]  = 

Coef  ficient [  fedu[  [  j  ]  ]  ,AAIi]  ] ,  {  j  ,Np''2> ,  {i,Np"2}  ] 

f edv=Expand [ edv ] 

Do[M[[j+Np"2,i]]  = 

Coefficient [f edu [  [  j]  ]  ,BB[i]  ] , {  j,Np''2} ,  {i,Np^2>] 
Do[M[[ j+Np"2,i+Np"2]]  = 

Coefficient [f edu [  [  j] ],CC[i] ] , {  j,Np^2}, {i,Np''2}] 
Do[M[  [  j+Np''2,2*Np''2]  ]  = 

Coefficient[fedu[  [  j]],AA[i]  ],{j,Np''2},{i,Np''2>] 

f edw=Expand [ edw ] 

Do[M[[2*Np"2+l,i]]  = 

Coefficient [f edu [  [  j]  ] ,BB[i]  ] , {  j,Np*2} ,  {i,Np''2} ] 
Do[M[[2*Np^2+l,i+Np^2]]  = 

Coef ficient (f edu [ [ j ] ] ,CC[i] ] , { j,Np^2> , {i,Np"2} ] 
DotM[ [2*Np^2+l,i+2*Np^2] ]  = 

Coef  ficient  [  f  edu  [  [  j  ]  ]  ,AA(i]  ] ,  {  j  ,Np*2} ,  {i,Np''2}  ] 

(*  Create  the  answer  vector  (bm)  for  the  m*sol=bm  system  *) 

bm  =  Array  [bbm, {2*Np''2+l,  1} ] 

Do  [bm[[i]]=0,  {i,2*Np''2}] 

Do  [bm[  [2*Np''2+l]  ]=Integrate  [P  ,{x,0,a},  {y,0,b}] 

(*  Solve  the  system  of  Equations  for  the  sol  vector  *) 
sol=LinearSolve [ m, bm ] 

(*  Determine  the  Bl,  Cl,  A1  coefficients  *) 

Do[BB[i]=sol[[i],{i,Np"2>] 

Do[CC[i]=sol[[i+Np''2},{i,Np"2}] 

Do[AA[i]=sol[  [i+2*Np^2]  ] ,  {i,Np''2}] 


(*  Creat  the  Displacement  functions  for  the  fiber  strata  *) 

Uf  [z_]=T2d>le[uff  [z]  [f  ] ,  {f  ,nplies}] 

Vf  [z_]=T£ible[vff  [z]  [f  ]  ,{f  ,nplies}] 

Wf [ z_]=Table[wf f [ z ] [ f ] , {f ,nplies>  j 

Do[uff  [z_]  [f  ]  =  Sum[zf  [  [f  ]  ]^i*phi[[i]  ],{i,Np''2},{f  ,nplies}] 


113 


Do[vff[z_][fl  =  Sum[z£[[f]ri*chi[[i]l,{i,Np^2},{f,nplies}] 
Do[wff  [z_]  [f  ]  =  Sum[  w[[l]  ],{i,Np''2},{f  ,nplies}] 

(*  Create  the  Displacement  Functions  for  the  matrix  strata  *) 

Um[z_]=Table[umm[z] [m] , {m,nplies+l} ] 
vm[z_]=Table[vmm[z] [m] , {m,nplies+l} ] 

Wm[ z_]=Table [ wmm[ z ] [mj , {m, nplies+1 } ] 

Do[umm[z_]  [f  ]  =  Sum[(al[[m]]+Bem[[m]]*z)''i*phi[[i]], 

{i,Np^2},{m, nplies+1}] 

Do[vmm[z_]  [f  ]  =  Sum[(al[[m]]+Bem[lm]]*z)''i*chi[[i]], 

{  i,Np'' 2  },  {m,  nplies+1 }  ] 

Do[wmm[z_] [f ]  =  Sum[  w[[lj], 

{i,Np^2},{m, nplies+1 } ] 

(*  Create  the  Strain  Functions  for  the  matrix  strata  *) 

efx[z_]=Table[effx[z] [f ] , {f ,nplies}] 
efy[z_]=Table[effy[z] [f ] , {f ,nplies} ] 
efz[z_]=Table[effz[z] [f ] , {f ,nplies}] 
efxy[z_]=Table[effxy[z] [f ] , {f ,nplies} ] 
efxz[z_]=Table[effxz[z] (f ] , {f ,nplies} ] 
efyz[z_]=Table[effyz[z] [f ] , {f ,nplies}] 

Do  [effx[z_] [f ]=D[Uf [z] [ [f 1 1 ,xl ,{f ,nplies}] 

DO  [effy[z_][f]=D[Vf[z][[f]],y],{f,nplies}] 

DO  [effz[2_][f]=0,{f,nplies>] 

DO  [effxlz_][f]=D[Uf[z][[f]],y]+D[Vf[z][[f]],x],{f,nplies}] 

DO  [effy[z_][f]*D[Uf[zl[[f]],z]+DtWf[2][tf]l,x],{f,nplies>] 

DO  [effy[z_][f]=D[Vf[z][[f]],zl+D[Wf[z][[f]3,y],{f,nplies>] 

(*  Create  the  Strain  Functions  for  the  matrix  strata  *) 

emx[z_]=Table[emmx[z] [m] , {m, nplies+1} ] 
emy[z_]=Table[emny(z] [m] , {m, nplies+1}] 
emz[z_]=Table[emmztz] [m] , {m, nplies+1} ] 
emxy[z_]=Table[emmxy[z] [m] , {m, nplies+1 } ] 
emxz[z_]=Table[emmxz[z] [m] , {m,nplies+l} ] 
emyz [ z_]=Table[emrayz [ z ] [m] , {m, nplies+1 } ] 

DO  [emmx[z_][f]=D[Um[z][[f]],x],{f,nplies}] 

Do  [emmy[z_][f]=D[Vm[z]IIf]],y],{f,nplies}] 

Do  [emmz[z_] [f ]=0, {f ,nplies} ] 

DO  [emmx[z_] [f ]=D[Um[z] [ [m] ] ,y l+D[Vf [z] [ [m] ] ,x] , {m, nplies+1} ] 

Do  [emn^[z_] [f ]=D[UmIz] [ [m] ] ,zl+D[Wf [z] [ [m] ] ,x] ,{m, nplies+1 } ] 

Do  [emny[z_] [f ]=D[Vm[z] [ [m] ] ,z]+D[Wf [z] [ [m] ] ,y] , {m, nplies+1}] 

(*  Create  the  Strain  vector  for  the  fiber  strata  *) 

Fstrain [  z_]  =Ted3le  [  f  strain [  z  ]  [  f  ] ,  { f ,  nplies }  ] 

Do[fstrain[z_] [f ]={efx[z] ,efy[z] [f ] ,efz[z] [f ] ,efxy[z] [f ] , 

efxz[z][f],efyz[z][f]},{f, nplies}] 

{*  Create  the  Strain  vector  for  the  matrix  strata  *) 
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Mstrain[z_]=Table[mstrain[z] [m] ,{m,nplies+l} ] 

Do[instrain[z_]  [f  ]={enix[z]  [m]  ,eny[z]  [m],einz[z]  [m]  ,einxy[z]  [m] , 

emxz [ z ] [ m] , emyz [ z ] [ m] } , {m, nplies+1 } ] 

(*  Calculate  the  Stress  in  each  stratum  by  Stiff ness*strain  *) 

Fstress[z_]=  Table[fstress[z] [f ],{f ,nplies}] 

Do[fstress[z_] [f ]=  S[ [f ] ] .Fstrain[z] [ [f ] ] , {f ,nplies}] 

Mstress[z_]=  Table[mstress[z] [m] , {m, nplies+1 } ] 

Do[mstress[z_] [m]=  S[ [m] ] .Fstrain[z] [ [m] ] , {m,nplies+l} ] 

(*  Place  the  stress  componentes  into  plotable  form  *) 

smx[z_]=Teible[smxx[z]  [m] ,  {m,nplies+l}  ] 
siiiy[z_]=Table[srayy[z]  [m] ,  {m,nplies+l>] 
smz[z_]=Tcible[smzz[z]  [m] ,  {m,nplies+l}  ] 
smxy [ z_] =Table [ smxyy [ z ] [ m ] , {m, nplies+1 } ] 
smxz [ z_] =Table [ smxzz [ z ] [ m ] , {m, nplies+1 } ] 
smyz[z_]=Table[smyzz[z] [m] » {m, nplies+1 } ] 

sfx[z_]=Table[sfxx[z] [f ] , {f ,nplies}] 
sfy[z_]=Tcible[sfyy[z]  [f  ] ,  {f  ,nplies}] 
sfz[z_]=Table[sfzz[z] [f ] , {f ,nplies}] 
sf xy ( z_] =Table [ sf xyy [ z ] [ f ] , {f , nplies} ] 
sfxz[z_]=Table[sfxzz[z] [f ], {f ,nplies}] 
sfyz[z_]=Table[sfyzz[z][f If /nplies}] 

Do[smxx[z_]  [m]  =sMstress [ z ]  [  [m]  ]  [  [1]  ],  {m, nplies+1 }  ] 

Dot  smyy [ z_] [m]  =Mstress [zi[[m]]tt2]], {m, nplies+1 } ] 

Dotsmzz[z_] [m]  =Mstress[z] [ [m] ] [ [3] ] , {m, nplies+1 } ] 

Do[ smxyy [z_] [m]=Mstress[z] [ [m] ] [ [4] ] r {m, nplies+1}] 

Do[smxzz[z_] [m]=Mstresstz] [ [m] ] [ [5] ] , {m,nplies+l} ] 

Do[snyzz[z_] [m]=Mstress[z] [ [m] ] [ [6] ] / {m,nplies+l} ] 

Dotsfxx[z_]  [f  ]  ==Fstress [ z ]  [  [f  ]  ]  [  [1]  ] ,  {f, nplies}  ] 

Do[sfyy[z_][f]  =Fstress[z][[f]] [[2] ],{f, nplies}] 

Do[sfzz[z_] [f ]  =Fstress[z] [ [f ] ] [ [3] ] , {f, nplies}] 

Do[sfxyy[z_] [f ]=Fstress[z] [ [f ] ] [ [4] ], {f/ nplies}] 

Dotsfxzz[z_] [f ]=Fstress[z] [[f]][t5]], {f, nplies} ] 

Do [ s  f yzz [ z_ ] [ f ] =Fstress [z][tf]J[[6]]/{f/ nplies } ] 
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